Module m_MEF90_LinAlg_class
   implicit none
   private
   public :: mef90Vect
   public :: mef90Mat 
   public :: mef90Tens4OS

!!!
!!!  
!!!  m_MEF90_LinAlg_class: Starting a OO implementation of the basic classes in m_MEF90_LinAlg
!!!
!!!  (c) 2020 Blaise Bourdin bourdin@lsu.edu
!!!

   type, abstract :: mef90Vect
   end type
   type, abstract :: mef90Mat
   end type
   type, abstract :: mef90Tens4OS
   end type
End module m_MEF90_LinAlg_class

Module m_MEF90_LinAlg
#include "finclude/petscdef.h"
   Use m_MEF90_LinAlg_class
   Use m_MEF90_Parameters
   Use m_MEF90_Utils
   Use petsc
   
   IMPLICIT NONE
 
   Type, Extends(mef90Vect) :: Vect2D
      PetscReal          :: X
      PetscReal          :: Y
   End Type Vect2D
 
   Type, Extends(mef90Vect) :: Vect3D
      PetscReal          :: X
      PetscReal          :: Y
      PetscReal          :: Z
   End Type Vect3D
 
   Type, Extends(mef90Mat) :: Mat2D
      PetscReal          :: XX
      PetscReal          :: XY
      PetscReal          :: YX
      PetscReal          :: YY
   End Type Mat2D
 
   Type, Extends(mef90Mat) :: MatS2D
      PetscReal          :: XX
      PetscReal          :: YY
      PetscReal          :: XY
   End Type MatS2D
 
   Type, Extends(mef90Mat) :: Mat3D
      PetscReal          :: XX
      PetscReal          :: XY
      PetscReal          :: XZ
      PetscReal          :: YX
      PetscReal          :: YY
      PetscReal          :: YZ
      PetscReal          :: ZX
      PetscReal          :: ZY
      PetscReal          :: ZZ
   End Type Mat3D
 
   Type, Extends(mef90Mat) :: MatS3D
      PetscReal          :: XX
      PetscReal          :: YY
      PetscReal          :: ZZ
      PetscReal          :: YZ
      PetscReal          :: XZ
      PetscReal          :: XY
   End Type MatS3D
 
 !! After much hesitation,
 !! - the terms are numbered in alphabetical order (i.e. XXYX and not XYXX) 
 !! - the terms are stored in alphabetical order
 !! 2014-07: Changed ordering to rows of the upper triangular part
 !!          and naming to be consistent with Voigt notations
   Type, Extends(mef90Tens4OS) :: Tens4OS2D
      PetscReal          :: XXXX,XXYY,XXXY
      PetscReal          ::      YYYY,YYXY
      PetscReal          ::           XYXY
      
   End Type Tens4OS2D
 
   Type, Extends(mef90Tens4OS) :: Tens4OS3D
      PetscReal          :: XXXX,XXYY,XXZZ,XXYZ,XXXZ,XXXY
      PetscReal          ::      YYYY,YYZZ,YYYZ,YYXZ,YYXY
      PetscReal          ::           ZZZZ,ZZYZ,ZZXZ,ZZXY
      PetscReal          ::                YZYZ,YZXZ,YZXY
      PetscReal          ::                     XZXZ,XZXY
      PetscReal          ::                          XYXY
   End Type Tens4OS3D
   
   Type(Vect2D),Parameter    :: MEF90Vect2De1 = Vect2D(1.0_Kr,0.0_Kr)
   Type(Vect2D),Parameter    :: MEF90Vect2De2 = Vect2D(0.0_Kr,1.0_Kr)

   Type(Vect3D),Parameter    :: MEF90Vect3De1 = Vect3D(1.0_Kr,0.0_Kr,0.0_Kr)
   Type(Vect3D),Parameter    :: MEF90Vect3De2 = Vect3D(0.0_Kr,1.0_Kr,0.0_Kr)
   Type(Vect3D),Parameter    :: MEF90Vect3De3 = Vect3D(0.0_Kr,0.0_Kr,1.0_Kr)
   
   Type(Mat2D),Parameter     :: MEF90Mat2DIdentity = Mat2D(1.0_Kr,0.0_Kr,          &
                                                           0.0_Kr,1.0_Kr)
   Type(MatS2D),Parameter    :: MEF90MatS2DIdentity = MatS2D(1.0_Kr,1.0_Kr,0.0_Kr)
   Type(Mat3D),Parameter     :: MEF90Mat3DIdentity = Mat3D(1.0_Kr,0.0_Kr,0.0_Kr,   &
                                                           0.0_Kr,1.0_Kr,0.0_Kr,   &
                                                           0.0_Kr,0.0_Kr,1.0_Kr)
   Type(MatS3D),Parameter    :: MEF90MatS3DIdentity = MatS3D(1.0_Kr,1.0_Kr,1.0_Kr, &
                                                             0.0_Kr,0.0_Kr,0.0_Kr)

   Type(Tens4OS2D),Parameter :: MEF90Tens4OS2DIdentity = Tens4OS2D(1.0_Kr,0.0_Kr,0.0_Kr, &
                                                                          1.0_Kr,0.0_Kr, & 
                                                                                 0.5_Kr)
   Type(Tens4OS3D),Parameter :: MEF90Tens4OS3DIdentity = Tens4OS3D(1.0_Kr,0.0_Kr,0.0_Kr,0.0_Kr,0.0_Kr,0.0_Kr, &
                                                                          1.0_Kr,0.0_Kr,0.0_Kr,0.0_Kr,0.0_Kr, &
                                                                                 1.0_Kr,0.0_Kr,0.0_Kr,0.0_Kr, &
                                                                                        0.5_Kr,0.0_Kr,0.0_Kr, &
                                                                                               0.5_Kr,0.0_Kr, &
                                                                                                      0.5_Kr)

   Interface Operator (+)
      Module Procedure SumVect2D,SumVect3D,SumMat2D,SumMat3D,SumMatS2D,SumMatS3D,SumTens4OS2D,SumTens4OS3D
   End Interface

   Interface Operator (-)
      Module Procedure DifVect2D,DifVect3D,DifMat2D,DifMat3D,DifMatS2D,DifMatS3D,DifTens4OS2D,DifTens4OS3D
   End Interface

   Interface Operator (*)
      Module Procedure DbleXVect2D,Vect2DXDble,DbleXVect3D,Vect3DXDble,                          &
         DbleXMat2D,Mat2DXDble,DbleXMat3D,Mat3DXDble,                                            &
         DbleXMatS2D,MatS2DXDble,DbleXMatS3D,MatS3DXDble,                                        &
         MatXVect2D,MatXVect3D,MatXVect2DS,MatXVect3DS,                                          &
         DbleXTens4OS2D,Tens4OS2DXDble,Tens4OS2DXMatS2D,Tens4OS2DXMat2D,                         &
         DbleXTens4OS3D,Tens4OS3DXDble,Tens4OS3DXMatS3D,Tens4OS3DXMat3D,                         &
         Mat2DXMat2D,Mat2DXMatS2D,MatS2DXMatS2D,Mat3DXMat3D,Mat3DXMatS3D,MatS3DXMatS3D,          &
         DotP2D,DotP3D
   End Interface

   Interface Operator (/)
      Module Procedure Vect2DQuot,Vect3DQuot,Mat2DQuot,Mat3DQuot,MatS2DQuot,MatS3DQuot,Tens4OS2DQuot,Tens4OS3DQuot
   End Interface

   Interface Operator (.DotP.)
      Module Procedure DotP2D,DotP3D,ContP2D,ContP3D,ContP2DS,ContP3DS,                         &
         Mat2DDotMatS2D,MatS2DDotMat2D,Mat3DDotMatS3D,MatS3DDotMat3D
   End Interface
   
   Interface Operator (.CrossP.)
      Module Procedure CrossP3D
   End Interface
   
   Interface Transpose
      Module Procedure Transpose2D,Transpose3D
   End Interface
   
   Interface Invert
      Module Procedure InvertMat2D,InvertMatS2D,InvertMat3D,InvertMatS3D,InvertTens4OS2D,InvertTens4OS3D
   End Interface

   Interface Operator (.TensP.)
      Module Procedure TensPVect2D,TensPVect3D
   End Interface
   
   Interface Operator (.SymP.)
      Module Procedure SymPVect2D,SymPVect3D,SymPMatS2D,SymPMatS3D
   End Interface
   
   Interface Operator (.oDot.)
      Module Procedure oDotMatS2D,oDotMatS3D
   End Interface
   
   Interface Trace
      Module Procedure Trace2D,Trace3D,Trace2DS,Trace3DS
   End Interface
   
   Interface Det
      Module Procedure DetMat2D,DetMatS2D,DetMat3D,DetMatS3D
   End Interface

   Interface Assignment (=)
      Module Procedure Vect2D_Get_Real,Vect3D_Get_Real,              &
         Vect2D_Get_VectR,Vect3D_Get_VectR,                          &
         Vect2DEQ,Vect3DEQ,Mat2D_Get_Real,Mat3D_Get_Real,            &
         Mat2DEQ,Mat3DEQ,MatS2D_Get_Real,MatS3D_Get_Real,            &
         MatS2DEQ,MatS3DEQ,MatS2D_Get_VectR,MatS3D_Get_VectR,        &
         VectR_Get_MatS2D,VectR_Get_MatS3D,                          &
         Tens4OS2D_Get_Real,Tens4OS3D_Get_Real,                      &
         Tens4OS2DToArray,ArrayToTens4OS2D,                          &
         Tens4OS3DToArray,ArrayToTens4OS3D,                          &
         Mat2DGetArray,Mat3DGetArray,                                &
         MatS2DGetArray,MatS3DGetArray,                              &
         ArrayGetMat2D,ArrayGetMat3D,                                &
         ArrayGetMatS2D,ArrayGetMatS3D,                              &
         MatS3DToMat3D,MatS2DToMat2D,                                &
         Mat3DToMatS3D,Mat2DToMatS2D,                                &
         Tens4OS2D2Array4,Tens4OS3D2Array4,                          &
         Array42Tens4OS2D,Array42Tens4OS3D
   End Interface
  
   Interface Symmetrize
      Module Procedure Symmetrize2D,Symmetrize3D
   End Interface
   
   Interface DeviatoricPart
      Module Procedure DeviatoricPart2D,DeviatoricPart2DS,DeviatoricPart3D,DeviatoricPart3DS
   End Interface

   Interface HydrostaticPart
      Module Procedure HydrostaticPart2D,HydrostaticPart2DS,HydrostaticPart3D,HydrostaticPart3DS
   End Interface

   Interface MatRaRt
      Module Procedure RaRtMat2D,RaRtMatS2D,RaRtMat3D,RaRtMatS3D
   End Interface

   Interface MatRtaR
      Module Procedure RtaRMat2D,RtaRMatS2D,RtaRMat3D,RtaRMatS3D
   End Interface

   Interface Norm
      Module Procedure Vect2DNorm,Vect3DNorm,Mat2DNorm,MatS2DNorm,Mat3DNorm,MatS3DNorm
      ! Tens4OS2DNorm,Tens4OS3DNorm
   End Interface
   
   Interface simplexNormal
      Module Procedure simplexNormal2D,simplexNormal3D
   End Interface simplexNormal
   
   Interface sqrt
      Module Procedure Tens4OS2DSquareRoot, Tens4OS3DSquareRoot
   End Interface
   
   Interface SpectralDecomposition
      Module Procedure MatS2DSpectralDecomposition,MatS3DSpectralDecomposition
   End Interface

   Interface Diagonalize
      !!! Diagonalize(A,P,D) returns P,D such that A = P D P^{-1}
      !!! The diagonal entries of D are sorted in increasing order.
      Module Procedure MatS2DEigenVectorValues,MatS3DEigenVectorValues
   End Interface

   Interface Moment
      Module Procedure Mat2DMoment, MatS2DMoment, Mat3DMoment, MatS3DMoment
   End Interface

   Interface Tens4OSTransform
      Module Procedure Tens4OS2DTransform, Tens4OS3DTransform
   End Interface

   !Interface MatSymToMat
   !   Module Procedure MatS3DToMat3D,MatS2DToMat2D
   !End Interface



!!$  Type(Vect2D),Parameter       :: e1_2D = (/ 1.0_Kr,0.0_Kr /)
!!$  Type(Vect2D),Parameter       :: e2_2D = (/ 0.0_Kr,1.0_Kr /)
!!$
!!$  Type(Vect3D),Parameter       :: e1_3D = (/ 1.0_Kr,0.0_Kr,0.0_Kr /)
!!$  Type(Vect3D),Parameter       :: e2_3D = (/ 0.0_Kr,1.0_Kr,0.0_Kr /)
!!$  Type(Vect3D),Parameter       :: e3_3D = (/ 0.0_Kr,0.0_Kr,1.0_Kr /)


Contains
   Function SumVect2D (V1,V2)
      Type(Vect2D),Intent(IN)                     :: V1
      Type(Vect2D),Intent(IN)                     :: V2
      Type(Vect2D)                                :: SumVect2D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      SumVect2D%X = V1%X + V2%X
      SumVect2D%Y = V1%Y + V2%Y
      flops = 2.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function SumVect2D

   Function SumVect3D (V1,V2)
      Type(Vect3D),Intent(IN)                     :: V1,V2
      Type(Vect3D)                                :: SumVect3D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      SumVect3D%X = V1%X + V2%X 
      SumVect3D%Y = V1%Y + V2%Y
      SumVect3D%Z = V1%Z + V2%Z
      flops = 3.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function SumVect3D

   Function SumMat2D (M1,M2)
      Type(Mat2D),Intent(IN)                      :: M1,M2
      Type(Mat2D)                                 :: SumMat2D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      SumMat2D%XX = M1%XX + M2%XX
      SumMat2D%XY = M1%XY + M2%XY
      SumMat2D%YX = M1%YX + M2%YX
      SumMat2D%YY = M1%YY + M2%YY
      flops = 4.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function SumMat2D
   
   Function SumMatS2D (M1,M2)
      Type(MatS2D),Intent(IN)                     :: M1,M2
      Type(MatS2D)                                :: SumMatS2D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      SumMatS2D%XX = M1%XX + M2%XX
      SumMatS2D%XY = M1%XY + M2%XY
      SumMatS2D%YY = M1%YY + M2%YY
      flops = 3.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function SumMatS2D
   
   Function SumMat3D (M1,M2)
      Type(Mat3D),Intent(IN)                      :: M1,M2
      Type(Mat3D)                                 :: SumMat3D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      SumMat3D%XX = M1%XX + M2%XX
      SumMat3D%XY = M1%XY + M2%XY
      SumMat3D%XZ = M1%XZ + M2%XZ
      SumMat3D%YX = M1%YX + M2%YX
      SumMat3D%YY = M1%YY + M2%YY
      SumMat3D%YZ = M1%YZ + M2%YZ
      SumMat3D%ZX = M1%ZX + M2%ZX
      SumMat3D%ZY = M1%ZY + M2%ZY
      SumMat3D%ZZ = M1%ZZ + M2%ZZ
      flops = 9.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function SumMat3D
   
   Function SumMatS3D (M1,M2)
      Type(MatS3D),Intent(IN)                     :: M1,M2
      Type(MatS3D)                                :: SumMatS3D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      SumMatS3D%XX = M1%XX + M2%XX
      SumMatS3D%YY = M1%YY + M2%YY
      SumMatS3D%ZZ = M1%ZZ + M2%ZZ
      SumMatS3D%YZ = M1%YZ + M2%YZ
      SumMatS3D%XZ = M1%XZ + M2%XZ
      SumMatS3D%XY = M1%XY + M2%XY
      flops = 6.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function SumMatS3D
   
   Function SumTens4OS2D (T1,T2)
      Type(Tens4OS2D),Intent(IN)                  :: T1,T2
      Type(Tens4OS2D)                             :: SumTens4OS2D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      SumTens4OS2D%XXXX = T1%XXXX + T2%XXXX
      SumTens4OS2D%XXYY = T1%XXYY + T2%XXYY
      SumTens4OS2D%XXXY = T1%XXXY + T2%XXXY
      
      SumTens4OS2D%YYYY = T1%YYYY + T2%YYYY
      SumTens4OS2D%YYXY = T1%YYXY + T2%YYXY

      SumTens4OS2D%XYXY = T1%XYXY + T2%XYXY
      
      flops = 6.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function SumTens4OS2D
   
   Function SumTens4OS3D (T1,T2)
      Type(Tens4OS3D),Intent(IN)                  :: T1,T2
      Type(Tens4OS3D)                             :: SumTens4OS3D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      SumTens4OS3D%XXXX = T1%XXXX + T2%XXXX  
      SumTens4OS3D%XXYY = T1%XXYY + T2%XXYY  
      SumTens4OS3D%XXZZ = T1%XXZZ + T2%XXZZ  
      SumTens4OS3D%XXYZ = T1%XXYZ + T2%XXYZ  
      SumTens4OS3D%XXXZ = T1%XXXZ + T2%XXXZ  
      SumTens4OS3D%XXXY = T1%XXXY + T2%XXXY
          
      SumTens4OS3D%YYYY = T1%YYYY + T2%YYYY  
      SumTens4OS3D%YYZZ = T1%YYZZ + T2%YYZZ  
      SumTens4OS3D%YYYZ = T1%YYYZ + T2%YYYZ  
      SumTens4OS3D%YYXZ = T1%YYXZ + T2%YYXZ  
      SumTens4OS3D%YYXY = T1%YYXY + T2%YYXY  
          
      SumTens4OS3D%ZZZZ = T1%ZZZZ + T2%ZZZZ  
      SumTens4OS3D%ZZYZ = T1%ZZYZ + T2%ZZYZ  
      SumTens4OS3D%ZZXZ = T1%ZZXZ + T2%ZZXZ  
      SumTens4OS3D%ZZXY = T1%ZZXY + T2%ZZXY  
          
      SumTens4OS3D%YZYZ = T1%YZYZ + T2%YZYZ 
      SumTens4OS3D%YZXZ = T1%YZXZ + T2%YZXZ  
      SumTens4OS3D%YZXY = T1%YZXY + T2%YZXY  
          
      SumTens4OS3D%XZXZ = T1%XZXZ + T2%XZXZ  
      SumTens4OS3D%XZXY = T1%XZXY + T2%XZXY  
          
      SumTens4OS3D%XYXY = T1%XYXY + T2%XYXY  
      flops = 21.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function SumTens4OS3D

   ! Overloading "-"
   Function DifVect2D (V1,V2)
      Type(Vect2D),Intent(IN)                     :: V1
      Type(Vect2D),Intent(IN)                     :: V2
      Type(Vect2D)                                :: DifVect2D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      DifVect2D%X = V1%X-V2%X
      DifVect2D%Y = V1%Y-V2%Y
      flops = 2.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function DifVect2D
   
   Function DifVect3D (V1,V2)
      Type(Vect3D),Intent(IN)                     :: V1,V2
      Type(Vect3D)                                :: DifVect3D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      DifVect3D%X = V1%X - V2%X 
      DifVect3D%Y = V1%Y - V2%Y
      DifVect3D%Z = V1%Z - V2%Z
      flops = 3.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function DifVect3D
   
   Function DifMat2D (M1,M2)
      Type(Mat2D),Intent(IN)                      :: M1,M2
      Type(Mat2D)                                 :: DifMat2D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      DifMat2D%XX = M1%XX - M2%XX
      DifMat2D%XY = M1%XY - M2%XY
      DifMat2D%YX = M1%YX - M2%YX
      DifMat2D%YY = M1%YY - M2%YY
      flops = 4.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function DifMat2D
   
   Function DifMatS2D (M1,M2)
      Type(MatS2D),Intent(IN)                     :: M1,M2
      Type(MatS2D)                                :: DifMatS2D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      DifMatS2D%XX = M1%XX - M2%XX
      DifMatS2D%YY = M1%YY - M2%YY
      DifMatS2D%XY = M1%XY - M2%XY
      flops = 3.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function DifMatS2D
   
   Function DifMat3D (M1,M2)
      Type(Mat3D),Intent(IN)                      :: M1,M2
      Type(Mat3D)                                 :: DifMat3D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      DifMat3D%XX = M1%XX - M2%XX
      DifMat3D%XY = M1%XY - M2%XY
      DifMat3D%XZ = M1%XZ - M2%XZ
      DifMat3D%YX = M1%YX - M2%YX
      DifMat3D%YY = M1%YY - M2%YY
      DifMat3D%YZ = M1%YZ - M2%YZ
      DifMat3D%ZX = M1%ZX - M2%ZX
      DifMat3D%ZY = M1%ZY - M2%ZY
      DifMat3D%ZZ = M1%ZZ - M2%ZZ
      flops = 9.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function DifMat3D
   
   Function DifMatS3D (M1,M2)
      Type(MatS3D),Intent(IN)                     :: M1,M2
      Type(MatS3D)                                :: DifMatS3D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      DifMatS3D%XX = M1%XX - M2%XX
      DifMatS3D%YY = M1%YY - M2%YY
      DifMatS3D%ZZ = M1%ZZ - M2%ZZ
      DifMatS3D%YZ = M1%YZ - M2%YZ
      DifMatS3D%XZ = M1%XZ - M2%XZ
      DifMatS3D%XY = M1%XY - M2%XY
      flops = 6.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function DifMatS3D
   
   Function DifTens4OS2D (T1,T2)
      Type(Tens4OS2D),Intent(IN)                  :: T1,T2
      Type(Tens4OS2D)                             :: DifTens4OS2D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      DifTens4OS2D%XXXX = T1%XXXX - T2%XXXX
      DifTens4OS2D%XXYY = T1%XXYY - T2%XXYY
      DifTens4OS2D%XXXY = T1%XXXY - T2%XXXY
      
      DifTens4OS2D%YYYY = T1%YYYY - T2%YYYY
      DifTens4OS2D%YYXY = T1%YYXY - T2%YYXY
      
      DifTens4OS2D%XYXY = T1%XYXY - T2%XYXY
      flops = 6.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function DifTens4OS2D
   
   Function DifTens4OS3D (T1,T2)
      Type(Tens4OS3D),Intent(IN)                  :: T1,T2
      Type(Tens4OS3D)                             :: DifTens4OS3D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      DifTens4OS3D%XXXX = T1%XXXX - T2%XXXX  
      DifTens4OS3D%XXYY = T1%XXYY - T2%XXYY  
      DifTens4OS3D%XXZZ = T1%XXZZ - T2%XXZZ  
      DifTens4OS3D%XXYZ = T1%XXYZ - T2%XXYZ  
      DifTens4OS3D%XXXZ = T1%XXXZ - T2%XXXZ  
      DifTens4OS3D%XXXY = T1%XXXY - T2%XXXY
      
      DifTens4OS3D%YYYY = T1%YYYY - T2%YYYY  
      DifTens4OS3D%YYZZ = T1%YYZZ - T2%YYZZ  
      DifTens4OS3D%YYYZ = T1%YYYZ - T2%YYYZ  
      DifTens4OS3D%YYXZ = T1%YYXZ - T2%YYXZ  
      DifTens4OS3D%YYXY = T1%YYXY - T2%YYXY  
      
      DifTens4OS3D%ZZZZ = T1%ZZZZ - T2%ZZZZ  
      DifTens4OS3D%ZZYZ = T1%ZZYZ - T2%ZZYZ  
      DifTens4OS3D%ZZXZ = T1%ZZXZ - T2%ZZXZ  
      DifTens4OS3D%ZZXY = T1%ZZXY - T2%ZZXY  
      
      DifTens4OS3D%YZYZ = T1%YZYZ - T2%YZYZ 
      DifTens4OS3D%YZXZ = T1%YZXZ - T2%YZXZ  
      DifTens4OS3D%YZXY = T1%YZXY - T2%YZXY  
      
      DifTens4OS3D%XZXZ = T1%XZXZ - T2%XZXZ  
      DifTens4OS3D%XZXY = T1%XZXY - T2%XZXY  
      
      DifTens4OS3D%XYXY = T1%XYXY - T2%XYXY  
      flops = 21.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function DifTens4OS3D

   ! Overloading "*"
   Function DbleXVect2D(D1,V1)
      PetscReal,Intent(IN)                        :: D1
      Type(Vect2D),Intent(IN)                     :: V1
      Type(Vect2D)                                :: DbleXVect2D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      DbleXVect2D%X = D1 * V1%X
      DbleXVect2D%Y = D1 * V1%Y
      flops = 2.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function DbleXVect2D
   
   Function Vect2DXDble(V1,D1)
      PetscReal,Intent(IN)                        :: D1
      Type(Vect2D),Intent(IN)                     :: V1
      Type(Vect2D)                                :: Vect2DXDble
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      Vect2DXDble%X = D1 * V1%X
      Vect2DXDble%Y = D1 * V1%Y
      flops = 2.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function Vect2DXDble
   
   Function DbleXVect3D(D1,V1)
      PetscReal,Intent(IN)                        :: D1
      Type(Vect3D),Intent(IN)                     :: V1
      Type(Vect3D)                                :: DbleXVect3D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      DbleXVect3D%X = D1 * V1%X
      DbleXVect3D%Y = D1 * V1%Y
      DbleXVect3D%Z = D1 * V1%Z
      flops = 3.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function DbleXVect3D
   
   Function Vect3DXDble(V1,D1)
      PetscReal,Intent(IN)                        :: D1
      Type(Vect3D),Intent(IN)                     :: V1
      Type(Vect3D)                                :: Vect3DXDble
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      Vect3DXDble%X = D1 * V1%X
      Vect3DXDble%Y = D1 * V1%Y
      Vect3DXDble%Z = D1 * V1%Z
      flops = 3.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function Vect3DXDble
   
   Function DbleXMat2D(D1,M1)
      PetscReal,Intent(IN)                        :: D1
      Type(Mat2D),Intent(IN)                      :: M1
      Type(Mat2D)                                 :: DbleXMat2D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      DbleXMat2D%XX = D1 * M1%XX
      DbleXMat2D%XY = D1 * M1%XY
      DbleXMat2D%YX = D1 * M1%YX
      DbleXMat2D%YY = D1 * M1%YY
      flops = 4.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function DbleXMat2D
   
   Function Mat2DXDble(M1,D1)
      PetscReal,Intent(IN)                        :: D1
      Type(Mat2D),Intent(IN)                      :: M1
      Type(Mat2D)                                 :: Mat2DXDble
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      Mat2DXDble%XX = D1 * M1%XX
      Mat2DXDble%XY = D1 * M1%XY
      Mat2DXDble%YX = D1 * M1%YX
      Mat2DXDble%YY = D1 * M1%YY
      flops = 4.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function Mat2DXDble
   
   Function DbleXMatS2D(D1,M1)
      PetscReal,Intent(IN)                        :: D1
      Type(MatS2D),Intent(IN)                     :: M1
      Type(MatS2D)                                :: DbleXMatS2D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      DbleXMatS2D%XX = D1 * M1%XX
      DbleXMatS2D%YY = D1 * M1%YY
      DbleXMatS2D%XY = D1 * M1%XY
      flops = 3.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function DbleXMatS2D
   
   Function MatS2DXDble(M1,D1)
      PetscReal,Intent(IN)                        :: D1
      Type(MatS2D),Intent(IN)                     :: M1
      Type(MatS2D)                                :: MatS2DXDble
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      MatS2DXDble%XX = D1 * M1%XX
      MatS2DXDble%YY = D1 * M1%YY
      MatS2DXDble%XY = D1 * M1%XY
      flops = 3.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function MatS2DXDble
   
   Function DbleXMat3D(D1,M1)
      PetscReal,Intent(IN)                        :: D1
      Type(Mat3D),Intent(IN)                      :: M1
      Type(Mat3D)                                 :: DbleXMat3D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      DbleXMat3D%XX = D1 * M1%XX
      DbleXMat3D%XY = D1 * M1%XY
      DbleXMat3D%XZ = D1 * M1%XZ
      DbleXMat3D%YX = D1 * M1%YX
      DbleXMat3D%YY = D1 * M1%YY
      DbleXMat3D%YZ = D1 * M1%YZ
      DbleXMat3D%ZX = D1 * M1%ZX
      DbleXMat3D%ZY = D1 * M1%ZY
      DbleXMat3D%ZZ = D1 * M1%ZZ
      flops = 9.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function DbleXMat3D
   
   Function Mat3DXDble(M1,D1)
      PetscReal,Intent(IN)                        :: D1
      Type(Mat3D),Intent(IN)                      :: M1
      Type(Mat3D)                                 :: Mat3DXDble
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      Mat3DXDble%XX = D1 * M1%XX
      Mat3DXDble%XY = D1 * M1%XY
      Mat3DXDble%XZ = D1 * M1%XZ
      Mat3DXDble%YX = D1 * M1%YX
      Mat3DXDble%YY = D1 * M1%YY
      Mat3DXDble%YZ = D1 * M1%YZ
      Mat3DXDble%ZX = D1 * M1%ZX
      Mat3DXDble%ZY = D1 * M1%ZY
      Mat3DXDble%ZZ = D1 * M1%ZZ
      flops = 9.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function Mat3DXDble
   
   Function DbleXMatS3D(D1,M1)
      PetscReal,Intent(IN)                        :: D1
      Type(MatS3D),Intent(IN)                     :: M1
      Type(MatS3D)                                :: DbleXMatS3D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      DbleXMatS3D%XX = D1 * M1%XX
      DbleXMatS3D%YY = D1 * M1%YY
      DbleXMatS3D%ZZ = D1 * M1%ZZ
      DbleXMatS3D%YZ = D1 * M1%YZ
      DbleXMatS3D%XZ = D1 * M1%XZ
      DbleXMatS3D%XY = D1 * M1%XY
      flops = 6.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function DbleXMatS3D
   
   Function MatS3DXDble(M1,D1)
      Real(Kind = Kr),Intent(IN)                  :: D1
      Type(MatS3D),Intent(IN)                     :: M1
      Type(MatS3D)                                :: MatS3DXDble
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      MatS3DXDble%XX = D1 * M1%XX
      MatS3DXDble%YY = D1 * M1%YY
      MatS3DXDble%ZZ = D1 * M1%ZZ
      MatS3DXDble%YZ = D1 * M1%YZ
      MatS3DXDble%XZ = D1 * M1%XZ
      MatS3DXDble%XY = D1 * M1%XY
      flops = 6.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function MatS3DXDble
   
   Function DbleXTens4OS2D (D1,T1)
      PetscReal,Intent(IN)                        :: D1
      Type(Tens4OS2D),Intent(IN)                  :: T1
      Type(Tens4OS2D)                             :: DbleXTens4OS2D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      DbleXTens4OS2D%XXXX = D1 * T1%XXXX
      DbleXTens4OS2D%XXXY = D1 * T1%XXXY
      DbleXTens4OS2D%XXYY = D1 * T1%XXYY
      
      DbleXTens4OS2D%XYXY = D1 * T1%XYXY
      DbleXTens4OS2D%YYXY = D1 * T1%YYXY
      
      DbleXTens4OS2D%YYYY = D1 * T1%YYYY
      flops = 6.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function DbleXTens4OS2D
   
   Function DbleXTens4OS3D (D1,T1)
      PetscReal,Intent(IN)                        :: D1
      Type(Tens4OS3D),Intent(IN)                  :: T1
      Type(Tens4OS3D)                             :: DbleXTens4OS3D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      DbleXTens4OS3D%XXXX = D1 * T1%XXXX  
      DbleXTens4OS3D%XXYY = D1 * T1%XXYY
      DbleXTens4OS3D%XXZZ = D1 * T1%XXZZ  
      DbleXTens4OS3D%XXYZ = D1 * T1%XXYZ  
      DbleXTens4OS3D%XXXZ = D1 * T1%XXXZ  
      DbleXTens4OS3D%XXXY = D1 * T1%XXXY  
          
      DbleXTens4OS3D%YYYY = D1 * T1%YYYY  
      DbleXTens4OS3D%YYZZ = D1 * T1%YYZZ  
      DbleXTens4OS3D%YYYZ = D1 * T1%YYYZ  
      DbleXTens4OS3D%YYXZ = D1 * T1%YYXZ  
      DbleXTens4OS3D%YYXY = D1 * T1%YYXY  
          
      DbleXTens4OS3D%ZZZZ = D1 * T1%ZZZZ  
      DbleXTens4OS3D%ZZYZ = D1 * T1%ZZYZ  
      DbleXTens4OS3D%ZZXZ = D1 * T1%ZZXZ  
      DbleXTens4OS3D%ZZXY = D1 * T1%ZZXY  
          
      DbleXTens4OS3D%YZYZ = D1 * T1%YZYZ  
      DbleXTens4OS3D%YZXZ = D1 * T1%YZXZ  
      DbleXTens4OS3D%YZXY = D1 * T1%YZXY  
          
      DbleXTens4OS3D%XZXZ = D1 * T1%XZXZ  
      DbleXTens4OS3D%XZXY = D1 * T1%XZXY  
          
      DbleXTens4OS3D%XYXY = D1 * T1%XYXY  
      flops = 21.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function DbleXTens4OS3D
   
   Function Tens4OS2DXDble (T1,D1)
      PetscReal,Intent(IN)                        :: D1
      Type(Tens4OS2D),Intent(IN)                  :: T1
      Type(Tens4OS2D)                             :: Tens4OS2DXDble
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      Tens4OS2DXDble%XXXX = D1 * T1%XXXX
      Tens4OS2DXDble%XXXY = D1 * T1%XXXY
      Tens4OS2DXDble%XXYY = D1 * T1%XXYY
      
      Tens4OS2DXDble%XYXY = D1 * T1%XYXY
      Tens4OS2DXDble%YYXY = D1 * T1%YYXY
      
      Tens4OS2DXDble%YYYY = D1 * T1%YYYY
      flops = 6.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function Tens4OS2DXDble
   
   Function Tens4OS3DXDble (T1,D1)
      Type(Tens4OS3D),Intent(IN)                  :: T1
      PetscReal,Intent(IN)                        :: D1
      Type(Tens4OS3D)                             :: Tens4OS3DXDble
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      Tens4OS3DXDble%XXXX = D1 * T1%XXXX  
      Tens4OS3DXDble%XXYY = D1 * T1%XXYY
      Tens4OS3DXDble%XXZZ = D1 * T1%XXZZ  
      Tens4OS3DXDble%XXYZ = D1 * T1%XXYZ  
      Tens4OS3DXDble%XXXZ = D1 * T1%XXXZ  
      Tens4OS3DXDble%XXXY = D1 * T1%XXXY  
          
      Tens4OS3DXDble%YYYY = D1 * T1%YYYY  
      Tens4OS3DXDble%YYZZ = D1 * T1%YYZZ  
      Tens4OS3DXDble%YYYZ = D1 * T1%YYYZ  
      Tens4OS3DXDble%YYXZ = D1 * T1%YYXZ  
      Tens4OS3DXDble%YYXY = D1 * T1%YYXY  
          
      Tens4OS3DXDble%ZZZZ = D1 * T1%ZZZZ  
      Tens4OS3DXDble%ZZYZ = D1 * T1%ZZYZ  
      Tens4OS3DXDble%ZZXZ = D1 * T1%ZZXZ  
      Tens4OS3DXDble%ZZXY = D1 * T1%ZZXY  
          
      Tens4OS3DXDble%YZYZ = D1 * T1%YZYZ  
      Tens4OS3DXDble%YZXZ = D1 * T1%YZXZ  
      Tens4OS3DXDble%YZXY = D1 * T1%YZXY  
          
      Tens4OS3DXDble%XZXZ = D1 * T1%XZXZ  
      Tens4OS3DXDble%XZXY = D1 * T1%XZXY  
          
      Tens4OS3DXDble%XYXY = D1 * T1%XYXY  
      flops = 21.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function Tens4OS3DXDble

   Function  MatXVect2D(M1,V1)
      Type(Mat2D),Intent(IN)                      :: M1
      Type(Vect2D),Intent(IN)                     :: V1
      Type(Vect2D)                                :: MatXVect2D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      MatXVect2D%X = M1%XX * V1%X + M1%XY * V1%Y
      MatXVect2D%Y = M1%YX * V1%X + M1%YY * V1%Y
      flops = 6.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function MatXVect2D
   
   Function  MatXVect2DS(M1,V1)
      Type(MatS2D),Intent(IN)                     :: M1
      Type(Vect2D),Intent(IN)                     :: V1
      Type(Vect2D)                               :: MatXVect2DS
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      MatXVect2DS%X = M1%XX * V1%X + M1%XY * V1%Y
      MatXVect2DS%Y = M1%XY * V1%X + M1%YY * V1%Y
      flops = 6.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function MatXVect2DS
   
   Function  MatXVect3D(M1,V1)
      Type(Mat3D),Intent(IN)                      :: M1
      Type(Vect3D),Intent(IN)                     :: V1
      Type(Vect3D)                                :: MatXVect3D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      MatXVect3D%X = M1%XX * V1%X + M1%XY * V1%Y + M1%XZ * V1%Z
      MatXVect3D%Y = M1%YX * V1%X + M1%YY * V1%Y + M1%YZ * V1%Z
      MatXVect3D%Z = M1%ZX * V1%X + M1%ZY * V1%Y + M1%ZZ * V1%Z
      flops = 15.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function MatXVect3D
   
   Function  MatXVect3DS(M1,V1)
      Type(MatS3D),Intent(IN)                     :: M1
      Type(Vect3D),Intent(IN)                     :: V1
      Type(Vect3D)                                :: MatXVect3DS
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      MatXVect3DS%X = M1%XX * V1%X + M1%XY * V1%Y + M1%XZ * V1%Z
      MatXVect3DS%Y = M1%XY * V1%X + M1%YY * V1%Y + M1%YZ * V1%Z
      MatXVect3DS%Z = M1%XZ * V1%X + M1%YZ * V1%Y + M1%ZZ * V1%Z
      flops = 15.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function MatXVect3DS
   
   Function Tens4OS2DXMatS2D(T1,M1)
      Type(Tens4OS2D),Intent(IN)                  :: T1
      Type(MatS2D),Intent(IN)                     :: M1
      Type(MatS2D)                                :: Tens4OS2DXMatS2D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      Tens4OS2DXMatS2D%XX = T1%XXXX * M1%XX + T1%XXYY * M1%YY + T1%XXXY * M1%XY * 2.0_Kr
      Tens4OS2DXMatS2D%YY = T1%XXYY * M1%XX + T1%YYYY * M1%YY + T1%YYXY * M1%XY * 2.0_Kr
      Tens4OS2DXMatS2D%XY = T1%XXXY * M1%XX + T1%YYXY * M1%YY + T1%XYXY * M1%XY * 2.0_Kr
      flops = 18.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function Tens4OS2DXMatS2D
   
   Function Tens4OS3DXMatS3D(T1,M1)
      Type(Tens4OS3D),Intent(IN)                  :: T1
      Type(MatS3D),Intent(IN)                     :: M1
      Type(MatS3D)                                :: Tens4OS3DXMatS3D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      Tens4OS3DXMatS3D%XX =  T1%XXXX * M1%XX + T1%XXYY * M1%YY + T1%XXZZ * M1%ZZ &
                          + (T1%XXYZ * M1%YZ + T1%XXXZ * M1%XZ + T1%XXXY * M1%XY) * 2.0_Kr
      Tens4OS3DXMatS3D%YY =  T1%XXYY * M1%XX + T1%YYYY * M1%YY + T1%YYZZ * M1%ZZ &
                          + (T1%YYYZ * M1%YZ + T1%YYXZ * M1%XZ + T1%YYXY * M1%XY) * 2.0_Kr
      Tens4OS3DXMatS3D%ZZ =  T1%XXZZ * M1%XX + T1%YYZZ * M1%YY + T1%ZZZZ * M1%ZZ &
                          + (T1%ZZYZ * M1%YZ + T1%ZZXZ * M1%XZ + T1%ZZXY * M1%XY) * 2.0_Kr
      Tens4OS3DXMatS3D%YZ =  T1%XXYZ * M1%XX + T1%YYYZ * M1%YY + T1%ZZYZ * M1%ZZ &
                          + (T1%YZYZ * M1%YZ + T1%YZXZ * M1%XZ + T1%YZXY * M1%XY) * 2.0_Kr
      Tens4OS3DXMatS3D%XZ =  T1%XXXZ * M1%XX + T1%YYXZ * M1%YY + T1%ZZXZ * M1%ZZ &
                          + (T1%YZXZ * M1%YZ + T1%XZXZ * M1%XZ + T1%XZXY * M1%XY) * 2.0_Kr
      Tens4OS3DXMatS3D%XY =  T1%XXXY * M1%XX + T1%YYXY * M1%YY + T1%ZZXY * M1%ZZ &
                          + (T1%YZXY * M1%YZ + T1%XZXY * M1%XZ + T1%XYXY * M1%XY) * 2.0_Kr
      flops = 66.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function Tens4OS3DXMatS3D
   
#undef __FUNCT__
#define __FUNCT__ "Tens4OS2DXMat2D"
!!!
!!!  
!!!  Tens4OS2DXMat2D:
!!!  
!!!  (c) 2016 Blaise Bourdin bourdin@lsu.edu
!!!
   Function Tens4OS2DXMat2D(T1,M1)
      Type(Tens4OS2D),Intent(IN)                  :: T1
      Type(Mat2D),Intent(IN)                      :: M1
      Type(MatS2D)                                :: Tens4OS2DXMat2D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      Tens4OS2DXMat2D = T1 * symmetrize(M1)
   End Function Tens4OS2DXMat2D
   
#undef __FUNCT__
#define __FUNCT__ "Tens4OS3DXMat3D"
!!!
!!!  
!!!  Tens4OS3DXMat3D:
!!!  
!!!  (c) 2016 Blaise Bourdin bourdin@lsu.edu
!!!
   Function Tens4OS3DXMat3D(T1,M1)
      Type(Tens4OS3D),Intent(IN)                  :: T1
      Type(Mat3D),Intent(IN)                      :: M1
      Type(MatS3D)                                :: Tens4OS3DXMat3D
      PetscInt                                    :: ierr
      
      Tens4OS3DXMat3D =  T1 * symmetrize(M1)
   End Function Tens4OS3DXMat3D
   
   Function Mat2DXMat2D(M1,M2)
      Type(Mat2D),Intent(IN)                      :: M1,M2
      Type(Mat2D)                                 :: Mat2DXMat2D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      Mat2DXMat2D%XX = M1%XX * M2%XX + M1%XY * M2%YX
      Mat2DXMat2D%XY = M1%XX * M2%XY + M1%XY * M2%YY
      Mat2DXMat2D%YX = M1%YX * M2%XX + M1%YY * M2%YX
      Mat2DXMat2D%YY = M1%YX * M2%XY + M1%YY * M2%YY
      flops = 12.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function Mat2DXMat2D

   Function Mat2DXMatS2D(M1,M2)
      Type(Mat2D),Intent(IN)                      :: M1
      Type(MatS2D),Intent(IN)                     :: M2
      Type(Mat2D)                                 :: Mat2DXMatS2D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      Mat2DXMatS2D%XX = M1%XX * M2%XX + M1%XY * M2%XY
      Mat2DXMatS2D%XY = M1%XX * M2%XY + M1%XY * M2%YY
      Mat2DXMatS2D%YX = M1%YX * M2%XX + M1%YY * M2%XY
      Mat2DXMatS2D%YY = M1%YX * M2%XY + M1%YY * M2%YY
      flops = 12.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function Mat2DXMatS2D
   
   Function MatS2DXMatS2D(M1,M2)
      Type(MatS2D),Intent(IN)                     :: M1,M2
      Type(Mat2D)                                 :: MatS2DXMatS2D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      MatS2DXMatS2D%XX = M1%XX * M2%XX + M1%XY * M2%XY
      MatS2DXMatS2D%XY = M1%XX * M2%XY + M1%XY * M2%YY
      MatS2DXMatS2D%YX = M1%XY * M2%XX + M1%YY * M2%XY
      MatS2DXMatS2D%YY = M1%XY * M2%XY + M1%YY * M2%YY
      flops = 12.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function MatS2DXMatS2D
      
   Function Mat3DXMat3D(M1,M2)
      Type(Mat3D),Intent(IN)                      :: M1,M2
      Type(Mat3D)                                 :: Mat3DXMat3D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      Mat3DXMat3D%XX = M1%XX * M2%XX + M1%XY * M2%YX + M1%XZ * M2%ZX
      Mat3DXMat3D%XY = M1%XX * M2%XY + M1%XY * M2%YY + M1%XZ * M2%ZY
      Mat3DXMat3D%XZ = M1%XX * M2%XZ + M1%XY * M2%YZ + M1%XZ * M2%ZZ 
      Mat3DXMat3D%YX = M1%YX * M2%XX + M1%YY * M2%YX + M1%YZ * M2%ZX
      Mat3DXMat3D%YY = M1%YX * M2%XY + M1%YY * M2%YY + M1%YZ * M2%ZY
      Mat3DXMat3D%YZ = M1%YX * M2%XZ + M1%YY * M2%YZ + M1%YZ * M2%ZZ
      Mat3DXMat3D%ZX = M1%ZX * M2%XX + M1%ZY * M2%YX + M1%ZZ * M2%ZX
      Mat3DXMat3D%ZY = M1%ZX * M2%XY + M1%ZY * M2%YY + M1%ZZ * M2%ZY
      Mat3DXMat3D%ZZ = M1%ZX * M2%XZ + M1%ZY * M2%YZ + M1%ZZ * M2%ZZ
      flops = 45.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function Mat3DXMat3D

   Function Mat3DXMatS3D(M1,M2)
      Type(Mat3D),Intent(IN)                      :: M1
      Type(MatS3D),Intent(IN)                     :: M2
      Type(Mat3D)                                 :: Mat3DXMatS3D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      Mat3DXMatS3D%XX = M1%XX * M2%XX + M1%XY * M2%XY + M1%XZ * M2%XZ
      Mat3DXMatS3D%XY = M1%XX * M2%XY + M1%XY * M2%YY + M1%XZ * M2%YZ
      Mat3DXMatS3D%XZ = M1%XX * M2%XZ + M1%XY * M2%YZ + M1%XZ * M2%ZZ 
      Mat3DXMatS3D%YX = M1%YX * M2%XX + M1%YY * M2%XY + M1%YZ * M2%XZ
      Mat3DXMatS3D%YY = M1%YX * M2%XY + M1%YY * M2%YY + M1%YZ * M2%YZ
      Mat3DXMatS3D%YZ = M1%YX * M2%XZ + M1%YY * M2%YZ + M1%YZ * M2%ZZ
      Mat3DXMatS3D%ZX = M1%ZX * M2%XX + M1%ZY * M2%XY + M1%ZZ * M2%XZ
      Mat3DXMatS3D%ZY = M1%ZX * M2%XY + M1%ZY * M2%YY + M1%ZZ * M2%YZ
      Mat3DXMatS3D%ZZ = M1%ZX * M2%XZ + M1%ZY * M2%YZ + M1%ZZ * M2%ZZ
      flops = 45.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function Mat3DXMatS3D
   
   Function MatS3DXMatS3D(M1,M2)
      Type(MatS3D),Intent(IN)                     :: M1,M2
      Type(Mat3D)                                 :: MatS3DXMatS3D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      MatS3DXMatS3D%XX = M1%XX * M2%XX + M1%XY * M2%XY + M1%XZ * M2%XZ
      MatS3DXMatS3D%XY = M1%XX * M2%XY + M1%XY * M2%YY + M1%XZ * M2%YZ
      MatS3DXMatS3D%XZ = M1%XX * M2%XZ + M1%XY * M2%YZ + M1%XZ * M2%ZZ 
      MatS3DXMatS3D%YX = M1%XY * M2%XX + M1%YY * M2%XY + M1%YZ * M2%XZ
      MatS3DXMatS3D%YY = M1%XY * M2%XY + M1%YY * M2%YY + M1%YZ * M2%YZ
      MatS3DXMatS3D%YZ = M1%XY * M2%XZ + M1%YY * M2%YZ + M1%YZ * M2%ZZ
      MatS3DXMatS3D%ZX = M1%XZ * M2%XX + M1%YZ * M2%XY + M1%ZZ * M2%XZ
      MatS3DXMatS3D%ZY = M1%XZ * M2%XY + M1%YZ * M2%YY + M1%ZZ * M2%YZ
      MatS3DXMatS3D%ZZ = M1%XZ * M2%XZ + M1%YZ * M2%YZ + M1%ZZ * M2%ZZ
      flops = 45.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function MatS3DXMatS3D
   
   ! Overloading "/"
   Function Vect2DQuot(V1,D1)
      Type(Vect2D),Intent(IN)                     :: V1
      PetscReal,Intent(IN)                        :: D1
      Type(Vect2D)                                :: Vect2DQuot
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      Vect2DQuot%X = V1%X / D1
      Vect2DQuot%Y = V1%Y / D1
      flops = 2.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function Vect2DQuot
   
   Function Vect3DQuot(V1,D1)
      Type(Vect3D),Intent(IN)                     :: V1
      PetscReal,Intent(IN)                        :: D1
      Type(Vect3D)                                :: Vect3DQuot
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      Vect3DQuot%X = V1%X / D1
      Vect3DQuot%Y = V1%Y / D1
      Vect3DQuot%Z = V1%Z / D1
      flops = 3.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function Vect3DQuot
   
   Function Mat2DQuot(M1,D1)
      PetscReal,Intent(IN)                        :: D1
      Type(Mat2D),Intent(IN)                      :: M1
      Type(Mat2D)                                 :: Mat2DQuot
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      Mat2DQuot%XX = M1%XX / D1
      Mat2DQuot%XY = M1%XY / D1
      Mat2DQuot%YX = M1%YX / D1
      Mat2DQuot%YY = M1%YY / D1
      flops = 4.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function Mat2DQuot
   
   Function MatS2DQuot(M1,D1)
      PetscReal,Intent(IN)                        :: D1
      Type(MatS2D),Intent(IN)                     :: M1
      Type(MatS2D)                                :: MatS2DQuot
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      MatS2DQuot%XX = M1%XX / D1
      MatS2DQuot%YY = M1%YY / D1
      MatS2DQuot%XY = M1%XY / D1
      flops = 3.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function MatS2DQuot
   
   Function Mat3DQuot(M1,D1)
      PetscReal,Intent(IN)                        :: D1
      Type(Mat3D),Intent(IN)                      :: M1
      Type(Mat3D)                                 :: Mat3DQuot
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      Mat3DQuot%XX = M1%XX / D1
      Mat3DQuot%XY = M1%XY / D1
      Mat3DQuot%XZ = M1%XZ / D1
      Mat3DQuot%YX = M1%YX / D1
      Mat3DQuot%YY = M1%YY / D1
      Mat3DQuot%YZ = M1%YZ / D1
      Mat3DQuot%ZX = M1%ZX / D1
      Mat3DQuot%ZY = M1%ZY / D1
      Mat3DQuot%ZZ = M1%ZZ / D1
      flops = 9.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function Mat3DQuot
   
   Function MatS3DQuot(M1,D1)
      PetscReal,Intent(IN)                        :: D1
      Type(MatS3D),Intent(IN)                     :: M1
      Type(MatS3D)                                :: MatS3DQuot
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      MatS3DQuot%XX = M1%XX / D1
      MatS3DQuot%YY = M1%YY / D1
      MatS3DQuot%ZZ = M1%ZZ / D1
      MatS3DQuot%YZ = M1%YZ / D1
      MatS3DQuot%XZ = M1%XZ / D1
      MatS3DQuot%XY = M1%XY / D1
      flops = 6.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function MatS3DQuot
   
   Function Tens4OS2DQuot(T1,D1)
      Type(Tens4OS2D),Intent(IN)                  :: T1
      PetscReal,Intent(IN)                        :: D1
      Type(Tens4OS2D)                             :: Tens4OS2DQuot
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      Tens4OS2DQuot%XXXX = T1%XXXX / D1
      Tens4OS2DQuot%XXXY = T1%XXXY / D1
      Tens4OS2DQuot%XXYY = T1%XXYY / D1
      
      Tens4OS2DQuot%XYXY = T1%XYXY / D1
      Tens4OS2DQuot%YYXY = T1%YYXY / D1
      
      Tens4OS2DQuot%YYYY = T1%YYYY / D1
      flops = 6.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function Tens4OS2DQuot
   
   Function Tens4OS3DQuot (T1,D1)
      Type(Tens4OS3D),Intent(IN)                  :: T1
      PetscReal,Intent(IN)                        :: D1
      Type(Tens4OS3D)                             :: Tens4OS3DQuot
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      Tens4OS3DQuot%XXXX = T1%XXXX / D1  
      Tens4OS3DQuot%XXYY = T1%XXYY / D1
      Tens4OS3DQuot%XXZZ = T1%XXZZ / D1  
      Tens4OS3DQuot%XXYZ = T1%XXYZ / D1  
      Tens4OS3DQuot%XXXZ = T1%XXXZ / D1  
      Tens4OS3DQuot%XXXY = T1%XXXY / D1  
          
      Tens4OS3DQuot%YYYY = T1%YYYY / D1  
      Tens4OS3DQuot%YYZZ = T1%YYZZ / D1  
      Tens4OS3DQuot%YYYZ = T1%YYYZ / D1  
      Tens4OS3DQuot%YYXZ = T1%YYXZ / D1  
      Tens4OS3DQuot%YYXY = T1%YYXY / D1  
          
      Tens4OS3DQuot%ZZZZ = T1%ZZZZ / D1  
      Tens4OS3DQuot%ZZYZ = T1%ZZYZ / D1  
      Tens4OS3DQuot%ZZXZ = T1%ZZXZ / D1  
      Tens4OS3DQuot%ZZXY = T1%ZZXY / D1  
          
      Tens4OS3DQuot%YZYZ = T1%YZYZ / D1  
      Tens4OS3DQuot%YZXZ = T1%YZXZ / D1  
      Tens4OS3DQuot%YZXY = T1%YZXY / D1  
          
      Tens4OS3DQuot%XZXZ = T1%XZXZ / D1  
      Tens4OS3DQuot%XZXY = T1%XZXY / D1  
          
      Tens4OS3DQuot%XYXY = T1%XYXY / D1  
      flops = 21.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function Tens4OS3DQuot
  
   ! dot product in 2D and 3D
   Function DotP2D(V1,V2)
      Type(Vect2D),Intent(IN)                     :: V1,V2
      PetscReal                                   :: DotP2D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      DotP2D = V1%X * V2%X + V1%Y * V2%Y
      flops = 3.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function DotP2D
   
   Function DotP3D(V1,V2)
      Type(Vect3D),Intent(IN)                     :: V1,V2
      PetscReal                                   :: DotP3D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      DotP3D = V1%X * V2%X + V1%Y * V2%Y + V1%Z * V2%Z
      flops = 5.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function DotP3D
   
   Function ContP2D(M1,M2)
      ! tr(A^t x B)
      Type(Mat2D),Intent(IN)                      :: M1,M2
      PetscReal                                   :: ContP2D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      ContP2D = M1%XX * M2%XX + M1%XY * M2%XY + M1%YX * M2%YX + M1%YY * M2%YY
      flops = 7.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function ContP2D
   
   Function ContP2DS(M1,M2)
      ! tr(A^t x B)
      Type(MatS2D),Intent(IN)                     :: M1,M2
      PetscReal                                   :: ContP2DS
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      ContP2DS = M1%XX * M2%XX + M1%YY * M2%YY + 2.0_Kr * M1%XY * M2%XY
      flops = 6.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function ContP2DS

#undef __FUNCT__
#define __FUNCT__ "Mat2DDotMatS2D"
!!!
!!!  
!!!  Mat2DDotMatS2D:
!!!  
!!!  (c) 2016 Blaise Bourdin bourdin@lsu.edu
!!!
   Function Mat2DDotMatS2D(M1,M2)
      Type(Mat2D),Intent(IN)                      :: M1
      Type(MatS2D),Intent(IN)                     :: M2
      PetscReal                                   :: Mat2DDotMatS2D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      Mat2DDotMatS2D = symmetrize(M1) .DotP. M2
   End Function Mat2DDotMatS2D

#undef __FUNCT__
#define __FUNCT__ "MatS2DDotMat2D"
!!!
!!!  
!!!  MatS2DDotMat2D:
!!!  
!!!  (c) 2016 Blaise Bourdin bourdin@lsu.edu
!!!
   Function MatS2DDotMat2D(M1,M2)
      Type(MatS2D),Intent(IN)                     :: M1
      Type(Mat2D),Intent(IN)                      :: M2
      PetscReal                                   :: MatS2DDotMat2D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      MatS2DDotMat2D = M1 .DotP. symmetrize(M2)
   End Function MatS2DDotMat2D

#undef __FUNCT__
#define __FUNCT__ "Mat2DDotMatS3D"
!!!
!!!  
!!!  Mat3DDotMatS3D:
!!!  
!!!  (c) 2016 Blaise Bourdin bourdin@lsu.edu
!!!
   Function Mat3DDotMatS3D(M1,M2)
      Type(Mat3D),Intent(IN)                      :: M1
      Type(MatS3D),Intent(IN)                     :: M2
      PetscReal                                   :: Mat3DDotMatS3D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      Mat3DDotMatS3D = symmetrize(M1) .DotP. M2
   End Function Mat3DDotMatS3D

#undef __FUNCT__
#define __FUNCT__ "MatS3DDotMat3D"
!!!
!!!  
!!!  MatS3DDotMat3D:
!!!  
!!!  (c) 2016 Blaise Bourdin bourdin@lsu.edu
!!!
   Function MatS3DDotMat3D(M1,M2)
      Type(MatS3D),Intent(IN)                     :: M1
      Type(Mat3D),Intent(IN)                      :: M2
      PetscReal                                   :: MatS3DDotMat3D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      MatS3DDotMat3D = M1 .DotP. symmetrize(M2)
   End Function MatS3DDotMat3D

   Function ContP3D(M1,M2)
      ! tr(A^t x B)
      Type(Mat3D),Intent(IN)                      :: M1,M2
      PetscReal                                   :: ContP3D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      ContP3D = M1%XX * M2%XX + M1%XY * M2%XY + M1%XZ * M2%XZ + &
           &    M1%YX * M2%YX + M1%YY * M2%YY + M1%YZ * M2%YZ + &
           &    M1%ZX * M2%ZX + M1%ZY * M2%ZY + M1%ZZ * M2%ZZ
      flops = 17.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function ContP3D
   
   Function ContP3DS(M1,M2)
      ! tr(A^t x B)
      Type(MatS3D),Intent(IN)                     :: M1,M2
      PetscReal                                   :: ContP3DS
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      ContP3DS = M1%XX * M2%XX + M1%YY * M2%YY + M1%ZZ * M2%ZZ + 2.0_Kr * M1%YZ * M2%YZ &
               + 2.0_Kr * M1%XZ * M2%XZ + 2.0_Kr * M1%XY * M2%XY
      flops = 14.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function ContP3DS
   
   ! cross product 3D
   Function CrossP3D(V1,V2)
      Type(Vect3D),Intent(IN)                     :: V1,V2
      Type(Vect3D)                                :: CrossP3D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      CrossP3D%X =  V1%Y * V2%Z - V1%Z * V2%Y
      CrossP3D%Y =  V1%Z * V2%X - V1%X * V2%Z
      CrossP3D%Z =  V1%X * V2%Y - V1%Y * V2%X
      flops = 9.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function CrossP3D
   
   ! Transpose
   Function Transpose2D(M1)
    Type(Mat2D),Intent(IN)                     :: M1
    Type(Mat2D)                                :: Transpose2D
   
    Transpose2D%XX = M1%XX
    Transpose2D%XY = M1%YX
    Transpose2D%YX = M1%XY
    Transpose2D%YY = M1%YY
   End Function Transpose2D
   
   Function Transpose3D(M1)
      Type(Mat3D),Intent(IN)                      :: M1
      Type(Mat3D)                                 :: Transpose3D
      
      Transpose3D%XX = M1%XX
      Transpose3D%XY = M1%YX
      Transpose3D%XZ = M1%ZX
      Transpose3D%YX = M1%XY
      Transpose3D%YY = M1%YY
      Transpose3D%YZ = M1%ZY
      Transpose3D%ZX = M1%XZ
      Transpose3D%ZY = M1%YZ
      Transpose3D%ZZ = M1%ZZ
   End Function Transpose3D

   
   ! Tensor product
   Function TensPVect2D (V1,V2)
      Type(Vect2D),Intent(IN)                     :: V1
      Type(Vect2D),Intent(IN)                     :: V2
      Type(Mat2D)                                 :: TensPVect2D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      TensPVect2D%XX = V1%X * V2%X
      TensPVect2D%XY = V1%X * V2%Y
      TensPVect2D%YX = V1%Y * V2%X
      TensPVect2D%YY = V1%Y * V2%Y
      flops = 4.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function TensPVect2D
   
   Function TensPVect3D (V1,V2)
      Type(Vect3D),Intent(IN)                     :: V1,V2
      Type(Mat3D)                                 :: TensPVect3D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      TensPVect3D%XX = V1%X * V2%X 
      TensPVect3D%XY = V1%X * V2%Y
      TensPVect3D%XZ = V1%X * V2%Z
      TensPVect3D%YX = V1%Y * V2%X 
      TensPVect3D%YY = V1%Y * V2%Y
      TensPVect3D%YZ = V1%Y * V2%Z
      TensPVect3D%ZX = V1%Z * V2%X 
      TensPVect3D%ZY = V1%Z * V2%Y
      TensPVect3D%ZZ = V1%Z * V2%Z
      flops = 9.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function TensPVect3D
   
   ! Symmetrized product
   Function SymPVect2D (V1,V2)
      Type(Vect2D),Intent(IN)                     :: V1,V2
      Type(MatS2D)                                :: SymPVect2D
      
      SymPVect2D = Symmetrize(V1 .TensP. V2)
   End Function SymPVect2D
   
   Function SymPVect3D (V1,V2)
      Type(Vect3D),Intent(IN)                     :: V1,V2
      Type(MatS3D)                                :: SymPVect3D
      
      SymPVect3D = Symmetrize(V1 .TensP. V2)
   End Function SymPVect3D
   
   Function SymPMatS2D(M1,M2)
      Type(MatS2D),Intent(IN)                     :: M1,M2
      Type(Tens4OS2D)                             :: SymPMatS2D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      SymPMatS2D%XXXX =  M1%XX * M2%XX
      SymPMatS2D%XXYY = (M1%XX * M2%YY + M1%YY * M2%XX) * 0.5_Kr
      SymPMatS2D%XXXY = (M1%XX * M2%XY + M1%XY * M2%XX) * 0.5_Kr

      SymPMatS2D%YYYY =  M1%YY * M2%YY
      SymPMatS2D%YYXY = (M1%YY * M2%XY + M1%XY * M2%YY) * 0.5_Kr

      SymPMatS2D%XYXY =  M1%XY * M2%XY
      flops = 15.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function SymPMatS2D
   
   Function SymPMatS3D(M1,M2)
      Type(MatS3D),Intent(IN)                     :: M1,M2
      Type(Tens4OS3D)                             :: SymPMatS3D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      SymPMatS3D%XXXX =  M1%XX * M2%XX
      SymPMatS3D%XXYY = (M1%XX * M2%YY + M1%YY * M2%XX) * 0.5_Kr
      SymPMatS3D%XXZZ = (M1%XX * M2%ZZ + M1%ZZ * M2%XX) * 0.5_Kr
      SymPMatS3D%XXYZ = (M1%XX * M2%YZ + M1%YZ * M2%XX) * 0.5_Kr
      SymPMatS3D%XXXZ = (M1%XX * M2%XZ + M1%XZ * M2%XX) * 0.5_Kr
      SymPMatS3D%XXXY = (M1%XX * M2%XY + M1%XY * M2%XX) * 0.5_Kr

      SymPMatS3D%YYYY =  M1%YY * M2%YY
      SymPMatS3D%YYZZ = (M1%YY * M2%ZZ + M1%ZZ * M2%YY) * 0.5_Kr
      SymPMatS3D%YYYZ = (M1%YY * M2%YZ + M1%YZ * M2%YY) * 0.5_Kr
      SymPMatS3D%YYXZ = (M1%YY * M2%XZ + M1%XZ * M2%YY) * 0.5_Kr
      SymPMatS3D%YYXY = (M1%YY * M2%XY + M1%XY * M2%YY) * 0.5_Kr

      SymPMatS3D%ZZZZ =  M1%ZZ * M2%ZZ
      SymPMatS3D%ZZYZ = (M1%ZZ * M2%YZ + M1%YZ * M2%ZZ) * 0.5_Kr
      SymPMatS3D%ZZXZ = (M1%ZZ * M2%XZ + M1%XZ * M2%ZZ) * 0.5_Kr
      SymPMatS3D%ZZXY = (M1%ZZ * M2%XY + M1%XY * M2%ZZ) * 0.5_Kr

      SymPMatS3D%YZYZ =  M1%YZ * M2%YZ
      SymPMatS3D%YZXZ = (M1%YZ * M2%XZ + M1%XZ * M2%YZ) * 0.5_Kr
      SymPMatS3D%YZXY = (M1%YZ * M2%XY + M1%XY * M2%YZ) * 0.5_Kr

      SymPMatS3D%XZXZ =  M1%XZ * M2%XZ
      SymPMatS3D%XZXY = (M1%XZ * M2%XY + M1%XY * M2%XZ) * 0.5_Kr

      SymPMatS3D%XYXY =  M1%XY * M2%XY
      flops = 66.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function SymPMatS3D


   Function oDotMatS2D(M1,M2)
      Type(MatS2D),Intent(IN)                     :: M1,M2
      Type(Tens4OS2D)                             :: oDotMatS2D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr

      oDotMatS2D%XXXX =  M1%XX + M2%XX
      oDotMatS2D%XXYY =  M1%XY * M2%XY
      oDotMatS2D%XXXY = (M1%XX * M2%XY + M1%XY * M2%XX) * 0.5_Kr

      oDotMatS2D%YYYY =  M1%YY + M2%YY
      oDotMatS2D%YYXY = (M1%XY * M2%YY + M1%YY * M2%XY) * 0.5_Kr
      
      oDotMatS2D%XYXY = (M1%XX * M2%YY + M1%YY * M2%XX) * 0.5_Kr
      flops = 15.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function oDotMatS2D


   Function oDotMatS3D(M1,M2)
      Type(MatS3D),Intent(IN)                     :: M1,M2
      Type(Tens4OS3D)                             :: oDotMatS3D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      oDotMatS3D%XXXX =  M1%XX * M2%XX
      oDotMatS3D%XXYY =  M1%XY * M2%XY
      oDotMatS3D%XXZZ =  M1%XZ * M2%XZ
      oDotMatS3D%XXYZ = (M1%XY * M2%XZ + M1%XZ * M2%YZ) * 0.5_Kr
      oDotMatS3D%XXXZ = (M1%XX * M2%XZ + M1%XZ * M2%XX) * 0.5_Kr
      oDotMatS3D%XXXY = (M1%XX * M2%XY + M1%XY * M2%XX) * 0.5_Kr

      oDotMatS3D%YYYY =  M1%YY * M2%YY
      oDotMatS3D%YYZZ =  M1%YZ * M2%YZ
      oDotMatS3D%YYYZ = (M1%YY * M2%YZ + M1%YZ * M2%YY) * 0.5_Kr
      oDotMatS3D%YYXZ = (M1%XY * M2%YZ + M1%YZ * M2%XY) * 0.5_Kr
      oDotMatS3D%YYXY = (M1%XY * M2%YY + M1%YY * M2%XY) * 0.5_Kr

      oDotMatS3D%ZZZZ =  M1%ZZ * M2%ZZ
      oDotMatS3D%ZZYZ = (M1%YZ * M2%ZZ + M1%ZZ * M2%YZ) * 0.5_Kr
      oDotMatS3D%ZZXZ = (M1%XZ * M2%ZZ + M1%ZZ * M2%XZ) * 0.5_Kr
      oDotMatS3D%ZZXY = (M1%XZ * M2%YZ + M1%YZ * M2%XZ) * 0.5_Kr

      oDotMatS3D%YZYZ = (M1%YY * M2%ZZ + M1%ZZ * M2%YY) * 0.5_Kr
      oDotMatS3D%YZXZ = (M1%XY * M2%ZZ + M1%ZZ * M2%XY) * 0.5_Kr
      oDotMatS3D%YZXY = (M1%XY * M2%YZ + M1%YZ * M2%XY) * 0.5_Kr

      oDotMatS3D%XZXZ = (M1%XX * M2%ZZ + M1%ZZ * M2%XX) * 0.5_Kr
      oDotMatS3D%XZXY = (M1%XX * M2%YZ + M1%YZ * M2%XX) * 0.5_Kr

      oDotMatS3D%XYXY = (M1%XX * M2%YY + M1%YY * M2%XX) * 0.5_Kr
      flops = 66.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function oDotMatS3D

   
   Function Trace2D(M1)
      Type(Mat2D),Intent(IN)                      :: M1
      PetscReal                                   :: Trace2D
      
      Trace2D = M1%XX + M1%YY
   End Function Trace2D
   
   Function Trace2DS(M1)
      Type(MatS2D),Intent(IN)                     :: M1
      PetscReal                                   :: Trace2DS
      
      Trace2DS = M1%XX + M1%YY
   End Function Trace2DS
   
   Function Trace3D(M1)
      Type(Mat3D),Intent(IN)                      :: M1
      PetscReal                                   :: Trace3D
      
      Trace3D = M1%XX + M1%YY + M1%ZZ
   End Function Trace3D
   
   Function Trace3DS(M1)
      Type(MatS3D),Intent(IN)                     :: M1
      PetscReal                                   :: Trace3DS
      
      Trace3DS = M1%XX + M1%YY + M1%ZZ
   End Function Trace3DS
   
   Subroutine Vect2D_Get_Real(V1,R1)
      Type(Vect2D),Intent(OUT)                    :: V1
      PetscReal,Intent(IN)                        :: R1
      
      V1%X = R1
      V1%Y = R1
   End Subroutine Vect2D_Get_Real
   
   Subroutine Vect3D_Get_Real(V1,R1)
      Type(Vect3D),Intent(OUT)                    :: V1
      PetscReal,Intent(IN)                        :: R1
      
      V1%X = R1
      V1%Y = R1
      V1%Z = R1
   End Subroutine Vect3D_Get_Real
   
   Subroutine Vect2D_Get_VectR(V1,R1)
      Type(Vect2D),Intent(OUT)                    :: V1
      PetscReal,      Dimension(2),Intent(IN)     :: R1
      
      V1%X = R1(1)
      V1%Y = R1(2)
   End Subroutine Vect2D_Get_VectR
   
   Subroutine Vect3D_Get_VectR(V1,R1)
      Type(Vect3D),Intent(OUT)                    :: V1
      PetscReal,      Dimension(3),Intent(IN)     :: R1
      
      V1%X = R1(1)
      V1%Y = R1(2)
      V1%Z = R1(3)
   End Subroutine Vect3D_Get_VectR
   
   Subroutine Vect2DEQ(V1,V2)
      Type(Vect2D),Intent(OUT)                    :: V1
      Type(Vect2D),Intent(IN)                     :: V2
      
      V1%X = V2%X
      V1%Y = V2%Y
   End Subroutine Vect2DEQ
   
   Subroutine Vect3DEQ(V1,V2)
      Type(Vect3D),Intent(OUT)                    :: V1
      Type(Vect3D),Intent(IN)                     :: V2
      
      V1%X = V2%X
      V1%Y = V2%Y
      V1%Z = V2%Z
   End Subroutine Vect3DEQ
   
   Subroutine Mat2D_Get_Real(M1,R1)
      Type(Mat2D),Intent(OUT)                     :: M1
      PetscReal,Intent(IN)                        :: R1
      
      M1%XX = R1 ;M1%XY = R1
      M1%YX = R1 ;M1%YY = R1
   End Subroutine Mat2D_Get_Real
   
   Subroutine Mat3D_Get_Real(M1,R1)
      Type(Mat3D),Intent(OUT)                     :: M1
      PetscReal,Intent(IN)                        :: R1
      
      M1%XX = R1 ;M1%XY = R1 ;M1%XZ = R1
      M1%YX = R1 ;M1%YY = R1 ;M1%YZ = R1
      M1%ZX = R1 ;M1%ZY = R1 ;M1%ZZ = R1
   End Subroutine Mat3D_Get_Real
   
   Subroutine Mat2DEQ(M1,M2)
      Type(Mat2D),Intent(OUT)                     :: M1
      Type(Mat2D),Intent(IN)                      :: M2
      
      M1%XX = M2%XX  ; M1%XY = M2%XY
      M1%YX = M2%YX  ; M1%YY = M2%YY
   End Subroutine Mat2DEQ
   
   Subroutine Mat3DEQ(M1,M2)
      Type(Mat3D),Intent(OUT)                     :: M1
      Type(Mat3D),Intent(IN)                      :: M2
      
      M1%XX = M2%XX  ; M1%XY = M2%XY  ; M1%XZ = M2%XZ
      M1%YX = M2%YX  ; M1%YY = M2%YY  ; M1%YZ = M2%YZ
      M1%ZX = M2%ZX  ; M1%ZY = M2%ZY  ; M1%ZZ = M2%ZZ
   End Subroutine Mat3DEQ
   
   Subroutine MatS2D_Get_Real(M1,R1)
      Type(MatS2D),Intent(OUT)                    :: M1
      PetscReal,Intent(IN)                        :: R1
      
      M1%XX = R1 ;M1%YY = R1 ;M1%XY = R1 
   End Subroutine MatS2D_Get_Real
   
   Subroutine MatS3D_Get_Real(M1,R1)
      Type(MatS3D),Intent(OUT)                    :: M1
      PetscReal,Intent(IN)                        :: R1
      
      M1%XX = R1 ;M1%YY = R1 ;M1%ZZ = R1
      M1%YZ = R1 ;M1%XZ = R1 ;M1%XY = R1
   End Subroutine MatS3D_Get_Real
   
   Subroutine MatS2DEQ(M1,M2)
      Type(MatS2D),Intent(OUT)                    :: M1
      Type(MatS2D),Intent(IN)                     :: M2
      
      M1%XX = M2%XX  ; M1%YY = M2%YY ;M1%XY = M2%XY
   End Subroutine MatS2DEQ
   
   Subroutine MatS3DEQ(M1,M2)
      Type(MatS3D),Intent(OUT)                    :: M1
      Type(MatS3D),Intent(IN)                     :: M2
      
      M1%XX = M2%XX  ; M1%YY = M2%YY  ; M1%ZZ = M2%ZZ
      M1%YZ = M2%YZ  ; M1%XZ = M2%XZ  ; M1%XY = M2%XY
   End Subroutine MatS3DEQ
   
   Subroutine MatS2D_Get_VectR(M1,R1)
      Type(MatS2D),Intent(OUT)                    :: M1
      PetscReal,Dimension(3),Intent(IN)           :: R1

      M1%XX = R1(1)
      M1%YY = R1(2)
      M1%XY = R1(3)
   End Subroutine MatS2D_Get_VectR
   
   Subroutine MatS3D_Get_VectR(M1,R1)
      Type(MatS3D),Intent(OUT)                    :: M1
      PetscReal,Dimension(6),Intent(IN)           :: R1

      M1%XX = R1(1)
      M1%YY = R1(2)
      M1%ZZ = R1(3)
      M1%YZ = R1(4)
      M1%XZ = R1(5)
      M1%XY = R1(6)
   End Subroutine MatS3D_Get_VectR
   
   Subroutine VectR_Get_MatS2D(R1,M1)
      PetscReal,Dimension(3),Intent(OUT)          :: R1
      Type(MatS2D),Intent(IN)                     :: M1

      R1(1) = M1%XX
      R1(2) = M1%YY
      R1(3) = M1%XY
   End Subroutine VectR_Get_MatS2D
   
   Subroutine VectR_Get_MatS3D(R1,M1)
      PetscReal,Dimension(6),Intent(OUT)          :: R1
      Type(MatS3D),Intent(IN)                     :: M1

      R1(1) = M1%XX
      R1(2) = M1%YY
      R1(3) = M1%ZZ
      R1(4) = M1%YZ
      R1(5) = M1%XZ
      R1(6) = M1%XY
   End Subroutine VectR_Get_MatS3D
   
   Subroutine Mat2DGetArray(M,A)
      Type(Mat2D),Intent(OUT)                     :: M
      PetscReal,Dimension(2,2),Intent(IN)         :: A
      
      M%XX = A(1,1)
      M%XY = A(1,2)
      M%YX = A(2,1)
      M%YY = A(2,2)
   End Subroutine Mat2DGetArray
   
   Subroutine MatS2DGetArray(M,A)
      Type(MatS2D),Intent(OUT)                    :: M
      PetscReal,Dimension(2,2),Intent(IN)         :: A
      
      M%XX = A(1,1)
      M%YY = A(2,2)
      M%XY = A(1,2)
   End Subroutine MatS2DGetArray
   
   Subroutine Mat3DGetArray(M,A)
      Type(Mat3D),Intent(OUT)                     :: M
      PetscReal,Dimension(3,3),Intent(IN)         :: A
      
      M%XX = A(1,1)
      M%XY = A(1,2)
      M%XZ = A(1,3)
      M%YX = A(2,1)
      M%YY = A(2,2)
      M%YZ = A(2,3)
      M%ZX = A(3,1)
      M%ZY = A(3,2)
      M%ZZ = A(3,3)
   End Subroutine Mat3DGetArray

   Subroutine MatS3DGetArray(M,A)
      Type(MatS3D),Intent(OUT)                    :: M
      PetscReal,Dimension(3,3),Intent(IN)         :: A
      
      M%XX = A(1,1)
      M%YY = A(2,2)
      M%ZZ = A(3,3)
      M%YZ = A(2,3)
      M%XZ = A(1,3)
      M%XY = A(1,2)
   End Subroutine MatS3DGetArray
   
   Subroutine ArrayGetMat2D(A,M)
      PetscReal,Dimension(2,2),Intent(OUT)        :: A
      Type(Mat2D),Intent(IN)                      :: M
      
      A(1,1) = M%XX 
      A(1,2) = M%XY 
      A(2,1) = M%YX 
      A(2,2) = M%YY 
   End Subroutine ArrayGetMat2D
   
   Subroutine ArrayGetMatS2D(A,M)
      PetscReal,Dimension(2,2),Intent(OUT)        :: A
      Type(MatS2D),Intent(IN)                     :: M
      
      A(1,1) = M%XX 
      A(1,2) = M%XY 
      A(2,1) = M%XY 
      A(2,2) = M%YY 
   End Subroutine ArrayGetMatS2D
   
   Subroutine ArrayGetMat3D(A,M)
      PetscReal,Dimension(3,3),Intent(OUT)        :: A
      Type(Mat3D),Intent(IN)                      :: M
      
      A(1,1) = M%XX
      A(1,2) = M%XY
      A(1,3) = M%XZ
      A(2,1) = M%YX
      A(2,2) = M%YY
      A(2,3) = M%YZ
      A(3,1) = M%ZX
      A(3,2) = M%ZY
      A(3,3) = M%ZZ
   End Subroutine ArrayGetMat3D
   
   Subroutine ArrayGetMatS3D(A,M)
      PetscReal,Dimension(3,3),Intent(OUT)        :: A
      Type(MatS3D),Intent(IN)                     :: M
      
      A(1,1) = M%XX
      A(1,2) = M%XY
      A(1,3) = M%XZ
      A(2,1) = M%XY
      A(2,2) = M%YY
      A(2,3) = M%YZ
      A(3,1) = M%XZ
      A(3,2) = M%YZ
      A(3,3) = M%ZZ
   End Subroutine ArrayGetMatS3D
   
   Subroutine Tens4OS2DEQ(T1,T2)
      Type(Tens4OS2D),Intent(OUT)                 :: T1
      Type(Tens4OS2D),Intent(IN)                  :: T2
      
      T1%XXXX = T2%XXXX
      T1%XXXY = T2%XXXY
      T1%XXYY = T2%XXYY
      
      T1%XYXY = T2%XYXY
      T1%YYXY = T2%YYXY
      
      T1%YYYY = T2%YYYY
   End Subroutine Tens4OS2DEQ
   
   Subroutine Tens4OS3DEQ(T1,T2)
      Type(Tens4OS3D),Intent(OUT)                 :: T1
      Type(Tens4OS3D),Intent(IN)                  :: T2
      
      T1%XXXX = T2%XXXX  
      T1%XXYY = T2%XXYY
      T1%XXZZ = T2%XXZZ  
      T1%XXYZ = T2%XXYZ  
      T1%XXXZ = T2%XXXZ  
      T1%XXXY = T2%XXXY  
         
      T1%YYYY = T2%YYYY  
      T1%YYZZ = T2%YYZZ
      T1%YYYZ = T2%YYYZ
      T1%YYXZ = T2%YYXZ
      T1%YYXY = T2%YYXY
         
      T1%ZZZZ = T2%ZZZZ  
      T1%ZZYZ = T2%ZZYZ
      T1%ZZXZ = T2%ZZXZ
      T1%ZZXY = T2%ZZXY
         
      T1%YZYZ = T2%YZYZ
      T1%YZXZ = T2%YZXZ
      T1%YZXY = T2%YZXY
         
      T1%XZXZ = T2%XZXZ  
      T1%XZXY = T2%XZXY
         
      T1%XYXY = T2%XYXY 
   End Subroutine Tens4OS3DEQ
   
   Subroutine Tens4OS2D_Get_Real(T1,D1)
      Type(Tens4OS2D),Intent(OUT)                 :: T1
      PetscReal,Intent(IN)                        :: D1
      
      T1%XXXX = D1
      T1%XXXY = D1
      T1%XXYY = D1
      
      T1%XYXY = D1
      T1%YYXY = D1
      
      T1%YYYY = D1
   End Subroutine Tens4OS2D_Get_Real
   
   Subroutine Tens4OS3D_Get_Real(T1,D1)
      Type(Tens4OS3D),Intent(OUT)                 :: T1
      PetscReal,Intent(IN)                        :: D1
      
      T1%XXXX = D1  
      T1%XXYY = D1
      T1%XXZZ = D1  
      T1%XXYZ = D1  
      T1%XXXZ = D1  
      T1%XXXY = D1  
         
      T1%YYYY = D1  
      T1%YYZZ = D1  
      T1%YYYZ = D1  
      T1%YYXZ = D1  
      T1%YYXY = D1  
         
      T1%ZZZZ = D1  
      T1%ZZYZ = D1  
      T1%ZZXZ = D1  
      T1%ZZXY = D1  
         
      T1%YZYZ = D1  
      T1%YZXZ = D1  
      T1%YZXY = D1  
         
      T1%XZXZ = D1  
      T1%XZXY = D1  
         
      T1%XYXY = D1  
   End Subroutine Tens4OS3D_Get_Real
  
   Subroutine Tens4OS2DToArray(A,T)
      PetscReal,Dimension(3,3),Intent(OUT)        :: A
      Type(Tens4OS2D),Intent(IN)                  :: T

      PetscInt                                    :: ierr
      
      A(1,1) = T%XXXX
      A(1,2) = T%XXYY
      A(1,3) = T%XXXY * 2.0_Kr

      A(2,1) = A(1,2)
      A(2,2) = T%YYYY
      A(2,3) = T%YYXY * 2.0_Kr

      A(3,1) = A(3,1)
      A(3,2) = A(2,3)
      A(3,3) = T%XYXY * 2.0_Kr
      Call PetscLogFlops(3.0_PFlop,ierr);CHKERRQ(ierr)
   End Subroutine Tens4OS2DToArray

   Subroutine ArrayToTens4OS2D(T,A)
      Type(Tens4OS2D),Intent(OUT)                 :: T
      PetscReal,Dimension(3,3),Intent(IN)         :: A

      PetscInt                                    :: ierr

      T%XXXX = A(1,1) 
      T%XXYY = A(1,2)
      T%XXXY = A(1,3) * .5_Kr
      
      T%YYYY = A(2,2)
      T%YYXY = A(2,3) * .5_Kr 

      T%XYXY = A(3,3) * .5_Kr 
      Call PetscLogFlops(3.0_PFlop,ierr);CHKERRQ(ierr)
   End Subroutine ArrayToTens4OS2D

   Subroutine Tens4OS3DToArray(A,T)
      PetscReal,Dimension(6,6),Intent(OUT)        :: A
      Type(Tens4OS3D),Intent(IN)                  :: T

      PetscInt                                    :: ierr
      
      A(1,1) = T%XXXX
      A(1,2) = T%XXYY
      A(1,3) = T%XXZZ
      A(1,4) = T%XXYZ * 2.0_Kr
      A(1,5) = T%XXXZ * 2.0_Kr
      A(1,6) = T%XXXY * 2.0_Kr

      A(2,1) = A(1,2)
      A(2,2) = T%YYYY
      A(2,3) = T%YYZZ
      A(2,4) = T%YYYZ * 2.0_Kr
      A(2,5) = T%YYXZ * 2.0_Kr
      A(2,6) = T%YYXY * 2.0_Kr

      A(3,1) = A(1,3)
      A(3,2) = A(2,3)
      A(3,3) = T%ZZZZ
      A(3,4) = T%ZZYZ * 2.0_Kr
      A(3,5) = T%ZZXZ * 2.0_Kr
      A(3,6) = T%ZZXY * 2.0_Kr
   
      A(4,1) = A(1,4)
      A(4,2) = A(2,4)
      A(4,3) = A(3,4)
      A(4,4) = T%YZYZ * 2.0_Kr
      A(4,5) = T%YZXZ * 2.0_Kr
      A(4,6) = T%YZXY * 2.0_Kr
   
      A(5,1) = A(1,5)
      A(5,2) = A(2,5)
      A(5,3) = A(3,5)
      A(5,4) = A(4,5)
      A(5,5) = T%XZXZ * 2.0_Kr
      A(5,6) = T%XZXY * 2.0_Kr
   
      A(6,1) = A(1,6)
      A(6,2) = A(2,6)
      A(6,3) = A(3,6)
      A(6,4) = A(4,6)
      A(6,5) = A(5,6)
      A(6,6) = T%XYXY * 2.0_Kr
      Call PetscLogFlops(15.0_PFlop,ierr);CHKERRQ(ierr)
   End Subroutine Tens4OS3DToArray

   Subroutine ArrayToTens4OS3D(T,A)
      Type(Tens4OS3D),Intent(OUT)                 :: T
      PetscReal,Dimension(6,6),Intent(IN)         :: A

      PetscInt                                    :: ierr

      T%XXXX = A(1,1)
      T%XXYY = A(1,2)
      T%XXZZ = A(1,3)
      T%XXYZ = A(1,4) * .5_Kr
      T%XXXZ = A(1,5) * .5_Kr
      T%XXXY = A(1,6) * .5_Kr

      T%YYYY = A(2,2) 
      T%YYZZ = A(2,3) 
      T%YYYZ = A(2,4) * .5_Kr
      T%YYXZ = A(2,5) * .5_Kr
      T%YYXY = A(2,6) * .5_Kr

      T%ZZZZ = A(3,3) 
      T%ZZYZ = A(3,4) * .5_Kr
      T%ZZXZ = A(3,5) * .5_Kr
      T%ZZXY = A(3,6) * .5_Kr
   
      T%YZYZ = A(4,4) * .5_Kr
      T%YZXZ = A(4,5) * .5_Kr
      T%YZXY = A(4,6) * .5_Kr
   
      T%XZXZ = A(5,5) * .5_Kr
      T%XZXY = A(5,6) * .5_Kr
   
      T%XYXY = A(6,6) * .5_Kr
      Call PetscLogFlops(15.0_PFlop,ierr);CHKERRQ(ierr)
   End Subroutine ArrayToTens4OS3D

   !!! Overloading euclidian norm of derived types
   PetscReal Function Vect2DNorm(V)
      Type(Vect2D),Intent(IN)                     :: V
      
      Vect2DNorm = sqrt(V%X**2 + V%Y**2)
   End Function Vect2DNorm
   
   PetscReal Function Vect3DNorm(V)
      Type(Vect3D),Intent(IN)                     :: V
      
      Vect3DNorm = sqrt(V%X**2 + V%Y**2 + V%Z**2)
   End Function Vect3DNorm
   
   PetscReal Function Mat2DNorm(M)
      Type(Mat2D),Intent(IN)                      :: M
      
      Mat2DNorm = sqrt(M%XX**2 + M%XY**2 + M%YX**2 + M%YY**2)
   End Function Mat2DNorm
   
   PetscReal Function MatS2DNorm(M)
      Type(MatS2D),Intent(IN)                     :: M
      
      MatS2DNorm = sqrt(M%XX**2 + 2.0_Kr * M%XY**2 + M%YY**2)
   End Function MatS2DNorm
   
   PetscReal Function Mat3DNorm(M)
      Type(Mat3D),Intent(IN)                      :: M
      
      Mat3DNorm = sqrt(M%XX**2 + M%XY**2 + M%XZ**2 + M%YX**2 + M%YY**2 + M%YZ**2 + M%ZX**2 + M%ZY**2 + M%ZZ**2)
   End Function Mat3DNorm
   
   PetscReal Function MatS3DNorm(M)
      Type(MatS3D),Intent(IN)                     :: M
      
      MatS3DNorm = sqrt(M%XX**2 + 2.0_Kr * M%XY**2 + 2.0_Kr * M%XZ**2 + M%YY**2 + 2.0_Kr * M%YZ**2 + M%ZZ**2 )
   End Function MatS3DNorm
   
   Function Symmetrize2D(M1)
      Type(Mat2D),Intent(IN)                      :: M1
      Type(MatS2D)                                :: Symmetrize2D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
                                      
      Symmetrize2D%XX = M1%XX
      Symmetrize2D%YY = M1%YY
      Symmetrize2D%XY = (M1%XY + M1%YX) * 0.5_Kr
      flops = 2.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function Symmetrize2D
   
   Function Symmetrize3D(M1)
      Type(Mat3D),Intent(IN)                      :: M1
      Type(MatS3D)                                :: Symmetrize3D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      Symmetrize3D%XX = M1%XX
      Symmetrize3D%YY = M1%YY
      Symmetrize3D%ZZ = M1%ZZ
      
      Symmetrize3D%YZ = (M1%YZ + M1%ZY) * 0.5_Kr
      Symmetrize3D%XZ = (M1%XZ + M1%ZX) * 0.5_Kr
      Symmetrize3D%XY = (M1%XY + M1%YX) * 0.5_Kr
      flops = 6.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function Symmetrize3D
   
   Subroutine MatS2DToMat2D(M1,M2)
      Type(Mat2D),Intent(OUT)                     :: M1
      Type(MatS2D),Intent(IN)                     :: M2

      M1%XX = M2%XX
      M1%XY = M2%XY
      M1%YX = M2%XY
      M1%YY = M2%YY
   End Subroutine MatS2DToMat2D

   Subroutine MatS3DToMat3D(M1,M2)
      Type(Mat3D),Intent(OUT)                     :: M1
      Type(MatS3D),Intent(IN)                     :: M2

      M1%XX = M2%XX
      M1%XY = M2%XY
      M1%XZ = M2%XZ
      M1%YX = M2%XY
      M1%YY = M2%YY
      M1%YZ = M2%YZ
      M1%ZX = M2%XZ
      M1%ZY = M2%YZ
      M1%ZZ = M2%ZZ
   End Subroutine MatS3DToMat3D

   Subroutine Mat2DToMatS2D(M1,M2)
      Type(MatS2D),Intent(OUT)                    :: M1
      Type(Mat2D),Intent(IN)                      :: M2

      M1%XX = M2%XX
      M1%XY = M2%XY
      M1%YY = M2%YY
   End Subroutine Mat2DToMatS2D

   Subroutine Mat3DToMatS3D(M1,M2)
      Type(MatS3D),Intent(OUT)                    :: M1
      Type(Mat3D),Intent(IN)                      :: M2

      M1%XX = M2%XX
      M1%YY = M2%YY
      M1%ZZ = M2%ZZ
      M1%YZ = M2%YZ
      M1%XZ = M2%XZ
      M1%XY = M2%XY
   End Subroutine Mat3DToMatS3D

   Function RARtMat2D(A,R)
      ! A <- R.A.R^T
      Type(Mat2D),Intent(IN)                      :: A
      Type(Mat2D),Intent(IN)                      :: R
      Type(Mat2D)                                 :: RARtMat2D

      RARtMat2D = R*A*transpose(R)
   End Function RaRtMat2D   

   Function RARtMatS2D(A,R)
      ! A <- R.A.R^T
      Type(MatS2D),Intent(IN)                     :: A
      Type(Mat2D),Intent(IN)                      :: R
      Type(MatS2D)                                :: RARtMatS2D

      RARtMatS2D = R*A*transpose(R)
   End Function RaRtMatS2D   

   Function RARtMat3D(A,R)
      ! A <- R.A.R^T
      Type(Mat3D),Intent(IN)                      :: A
      Type(Mat3D),Intent(IN)                      :: R
      Type(Mat3D)                                 :: RARtMat3D

      RARtMat3D = R*A*transpose(R)
   End Function RaRtMat3D   

   Function RARtMatS3D(A,R)
      ! A <- R.A.R^T
      Type(MatS3D),Intent(IN)                     :: A
      Type(Mat3D),Intent(IN)                      :: R
      Type(MatS3D)                                :: RARtMatS3D

      RARtMatS3D = R*A*transpose(R)
   End Function RaRtMatS3D   

   Function RtARMat2D(A,R)
      ! A <- Rt.A.R
      Type(Mat2D),Intent(IN)                      :: A
      Type(Mat2D),Intent(IN)                      :: R
      Type(Mat2D)                                 :: RtARMat2D

      RtARMat2D = transpose(R)*A*R
   End Function RtaRMat2D   

   Function RtaRMatS2D(A,R)
      ! A <- R^T.A.R
      Type(MatS2D),Intent(IN)                     :: A
      Type(Mat2D),Intent(IN)                      :: R
      Type(MatS2D)                                :: RtaRMatS2D

      RtaRMatS2D = transpose(R)*A*R
   End Function RtaRMatS2D   

   Function RtaRMat3D(A,R)
      ! A <- R^T.A.R
      Type(Mat3D),Intent(IN)                      :: A
      Type(Mat3D),Intent(IN)                      :: R
      Type(Mat3D)                                 :: RtaRMat3D

      RtaRMat3D = transpose(R)*A*R
   End Function RtaRMat3D   

   Function RtaRMatS3D(A,R)
      ! A <- R^T.A.R
      Type(MatS3D),Intent(IN)                     :: A
      Type(Mat3D),Intent(IN)                      :: R
      Type(MatS3D)                                :: RtaRMatS3D

      RtaRMatS3D = transpose(R)*A*R
   End Function RtaRMatS3D   

   Function DeviatoricPart2D(M1)
      Type(Mat2D),Intent(IN)                      :: M1
      Type(Mat2D)                                 :: DeviatoricPart2D
      
      PetscReal                                   :: M1_Trace
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      M1_Trace = Trace(M1)
      DeviatoricPart2D%XX = M1%XX - M1_Trace * 0.5_Kr
      DeviatoricPart2D%XY = M1%XY
      DeviatoricPart2D%YX = M1%YX
      DeviatoricPart2D%YY = M1%YY - M1_Trace * 0.5_Kr
      flops = 4.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function DeviatoricPart2D


   Function DeviatoricPart2DS(M1)
      Type(MatS2D),Intent(IN)                     :: M1
      Type(MatS2D)                                :: DeviatoricPart2DS
      
      PetscReal                                   :: M1_Trace
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      M1_Trace = Trace(M1)
      DeviatoricPart2DS%XX = M1%XX - M1_Trace * 0.5_Kr
      DeviatoricPart2DS%YY = M1%YY - M1_Trace * 0.5_Kr
      DeviatoricPart2DS%XY = M1%XY
      flops = 4.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function DeviatoricPart2DS

   Function DeviatoricPart3D(M1)
      Type(Mat3D),Intent(IN)                      :: M1
      Type(Mat3D)                                 :: DeviatoricPart3D
      
      PetscReal                                   :: M1_Trace
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      M1_Trace = Trace(M1)
      DeviatoricPart3D%XX = M1%XX - M1_Trace / 3.0_Kr
      DeviatoricPart3D%XY = M1%XY
      DeviatoricPart3D%XZ = M1%XZ

      DeviatoricPart3D%YX = M1%YX
      DeviatoricPart3D%YY = M1%YY - M1_Trace / 3.0_Kr
      DeviatoricPart3D%YZ = M1%YZ

      DeviatoricPart3D%ZX = M1%ZX
      DeviatoricPart3D%ZY = M1%ZY
      DeviatoricPart3D%ZZ = M1%ZZ - M1_Trace / 3.0_Kr
      flops = 6.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function DeviatoricPart3D


   Function DeviatoricPart3DS(M1)
      Type(MatS3D),Intent(IN)                     :: M1
      Type(MatS3D)                                :: DeviatoricPart3DS
      
      PetscReal                                   :: M1_Trace
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      M1_Trace = Trace(M1)
      DeviatoricPart3DS%XX = M1%XX - M1_Trace / 3.0_Kr
      DeviatoricPart3DS%YY = M1%YY - M1_Trace / 3.0_Kr
      DeviatoricPart3DS%ZZ = M1%ZZ - M1_Trace / 3.0_Kr

      DeviatoricPart3DS%YZ = M1%YZ
      DeviatoricPart3DS%XZ = M1%XZ
      DeviatoricPart3DS%XY = M1%XY
      flops = 6.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function DeviatoricPart3DS

   Function HydrostaticPart2D(M1)
      Type(Mat2D),Intent(IN)                      :: M1
      Type(Mat2D)                                 :: HydrostaticPart2D
      
      PetscReal                                   :: M1_Trace
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      M1_Trace = Trace(M1)
      HydrostaticPart2D%XX = M1_Trace * 0.5_Kr
      HydrostaticPart2D%XY = 0.0_Kr
      HydrostaticPart2D%YX = 0.0_Kr
      HydrostaticPart2D%YY = M1_Trace * 0.5_Kr
      flops = 2.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function HydrostaticPart2D


   Function HydrostaticPart2DS(M1)
      Type(MatS2D),Intent(IN)                     :: M1
      Type(MatS2D)                                :: HydrostaticPart2DS
      
      PetscReal                                   :: M1_Trace
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      M1_Trace = Trace(M1)
      HydrostaticPart2DS%XX = M1_Trace * 0.5_Kr
      HydrostaticPart2DS%YY = M1_Trace * 0.5_Kr
      HydrostaticPart2DS%XY = 0.0_Kr
      flops = 2.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function HydrostaticPart2DS

   Function HydrostaticPart3D(M1)
      Type(Mat3D),Intent(IN)                      :: M1
      Type(Mat3D)                                 :: HydrostaticPart3D
      
      PetscReal                                   :: M1_Trace
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      M1_Trace = Trace(M1)
      HydrostaticPart3D%XX = M1_Trace / 3.0_Kr
      HydrostaticPart3D%XY = 0.0_Kr
      HydrostaticPart3D%XZ = 0.0_Kr

      HydrostaticPart3D%YX = 0.0_Kr
      HydrostaticPart3D%YY = M1_Trace / 3.0_Kr
      HydrostaticPart3D%YZ = 0.0_Kr
                         
      HydrostaticPart3D%ZX = 0.0_Kr
      HydrostaticPart3D%ZY = 0.0_Kr
      HydrostaticPart3D%ZZ = M1_Trace / 3.0_Kr
      flops = 6.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function HydrostaticPart3D


   Function HydrostaticPart3DS(M1)
      Type(MatS3D),Intent(IN)                     :: M1
      Type(MatS3D)                                :: HydrostaticPart3DS
      
      PetscReal                                   :: M1_Trace
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      M1_Trace = Trace(M1)
      HydrostaticPart3DS%XX = M1_Trace / 3.0_Kr
      HydrostaticPart3DS%YY = M1_Trace / 3.0_Kr
      HydrostaticPart3DS%ZZ = M1_Trace / 3.0_Kr

      HydrostaticPart3DS%YZ = 0.0_Kr
      HydrostaticPart3DS%XZ = 0.0_Kr
      HydrostaticPart3DS%XY = 0.0_Kr
      flops = 6.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function HydrostaticPart3DS

!====================================================================
!             END OF OPERATOR OVERLOADING
!====================================================================

   Function Vol_Tetra_3D(V1,V2,V3,V4)
      Type(Vect3D),Intent(IN)                     :: V1,V2,V3,V4
      PetscReal                                   :: Vol_Tetra_3D
      
      Type(Vect3D)                                :: C1,C2,C3
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      C1 = V1 - V4
      C2 = V2 - V4
      C3 = V3 - V4
      Vol_Tetra_3D = ABS(C1%X * (C2%Y * C3%Z - C2%Z * C3%Y) - C1%Y * (C2%X * C3%Z - C2%Z * C3%X)  &
                       + C1%Z * (C2%X * C3%Y - C2%Y * C3%X) ) / 6.0_Kr
      flops = 18.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function Vol_Tetra_3D
   
   Function Area_Tri_2D(S1,S2,S3)
      Type(Vect2D),Intent(IN)                     :: S1,S2,S3
      PetscReal                                   :: Area_Tri_2D
      
      Type(Vect2D)                                :: C1,C2
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      C1 = S2-S1
      C2 = S3-S1
      Area_Tri_2D = ABS(C1%X * C2%Y - C1%Y * C2%X) * 0.5_Kr
      flops = 4.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function Area_Tri_2D
   
   Function Ht_Min_Tri_2D(S1,S2,S3)
      Type(Vect2D),Intent(IN)                     :: S1,S2,S3
      PetscReal                                   :: Ht_Min_Tri_2D
      
      Type(Vect2D)                                :: C1,C2,C3
      PetscReal                                   :: H1,H2,H3,AreaX2
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      C1 = S2-S3
      C2 = S3-S1
      C3 = S1-S2
      
      AreaX2 = abs(C1%X * C2%Y - C1%Y * C2%X)
      
      H1 = AreaX2 / SQRT( (C1 .DotP. C1) )
      H2 = AreaX2 / SQRT( (C2 .DotP. C2) )
      H3 = AreaX2 / SQRT( (C3 .DotP. C3) )
      
      Ht_Min_Tri_2D = Min (H1,H2,H3)
      flops = 12.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function Ht_Min_Tri_2D
   
   Function DetMat2D(M)
      Type(Mat2D),Intent(IN)                      :: M
      PetscReal                                   :: DetMat2D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      DetMat2D = M%XX * M%YY - M%XY * M%YX
      flops = 3.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function DetMat2D
   
   Function DetMatS2D(M)
      Type(MatS2D),Intent(IN)                     :: M
      PetscReal                                   :: DetMatS2D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      DetMatS2D = M%XX * M%YY - M%XY * M%XY
      flops = 3.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function DetMatS2D
   
   
   Function DetMat3D(M)
      Type(Mat3D),Intent(IN)                      :: M
      PetscReal                                   :: DetMat3D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      DetMat3D = M%XX * (M%YY * M%ZZ - M%ZY * M%YZ) - M%YX * (M%XY * M%ZZ - M%ZY * M%XZ) + M%ZX * (M%XY * M%YZ - M%YY * M%XZ)
      flops = 14.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function DetMat3D
   
   Function DetMatS3D(M)
      Type(MatS3D),Intent(IN)                     :: M
      PetscReal                                   :: DetMatS3D
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      DetMatS3D = M%XX * (M%YY * M%ZZ - M%YZ * M%YZ) - M%XY * (M%XY * M%ZZ - M%YZ * M%XZ) + M%XZ * (M%XY * M%YZ - M%YY * M%XZ)
      flops = 14.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function DetMatS3D

   Function InvertMat2D(M)
      Type(Mat2D),Intent(IN)                      :: M
      Type(Mat2D)                                 :: InvertMat2D
      
      Type(Mat2D)                                 :: CofMt
      PetscReal                                   :: DetM
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      DetM = M%XX * M%YY - M%XY * M%YX
      CofMt%XX =  M%YY
      CofMt%XY = -M%XY
      CofMt%YX = -M%YX
      CofMt%YY =  M%XX
      
      InvertMat2D = CofMt / DetM
      flops = 6.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function InvertMat2D

   Function InvertMatS2D(M)
      Type(MatS2D),Intent(IN)                     :: M
      Type(MatS2D)                                :: InvertMatS2D
      
      Type(MatS2D)                                :: CofMt
      PetscReal                                   :: DetM
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      DetM = M%XX * M%YY - M%XY **2
      CofMt%XX =  M%YY
      CofMt%XY = -M%XY
      CofMt%YY =  M%XX
      
      InvertMatS2D = CofMt / DetM
      flops = 5.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function InvertMatS2D

   Function InvertMat3D(M)
      Type(Mat3D),Intent(IN)                      :: M
      Type(Mat3D)                                 :: InvertMat3D
      
      Type(Mat3D)                                 :: CofMt
      PetscReal                                   :: DetM
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      DetM = M%XX * (M%YY * M%ZZ - M%ZY * M%YZ) - M%YX * (M%XY * M%ZZ - M%ZY * M%XZ) + M%ZX * (M%XY * M%YZ - M%YY * M%XZ)
           
      CofMt%XX =   M%YY * M%ZZ - M%ZY * M%YZ
      CofMt%YX = -(M%YX * M%ZZ - M%ZX * M%YZ)
      CofMt%ZX =   M%YX * M%ZY - M%ZX * M%YY
      CofMt%XY = -(M%XY * M%ZZ - M%ZY * M%XZ)
      CofMt%YY =   M%XX * M%ZZ - M%ZX * M%XZ
      CofMt%ZY = -(M%XX * M%ZY - M%ZX * M%XY)
      CofMt%XZ =   M%XY * M%YZ - M%YY * M%XZ
      CofMt%YZ = -(M%XX * M%YZ - M%YX * M%XZ)
      CofMt%ZZ =   M%XX * M%YY - M%YX * M%XY
       
      InvertMat3D = CofMt / DetM
      flops = 46.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function InvertMat3D
   
   Function InvertMatS3D(M)
      Type(MatS3D),Intent(IN)                     :: M
      Type(MatS3D)                                :: InvertMatS3D
      
      Type(MatS3D)                                :: CofMt
      PetscReal                                   :: DetM
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      
      DetM = M%XX * (M%YY * M%ZZ - M%YZ * M%YZ) - M%XY * (M%XY * M%ZZ - M%YZ * M%XZ) + M%XZ * (M%XY * M%YZ - M%YY * M%XZ)
           
      CofMt%XX =   M%YY * M%ZZ - M%YZ * M%YZ
      CofMt%XY = -(M%XY * M%ZZ - M%XZ * M%YZ)
      CofMt%XZ =   M%XY * M%YZ - M%XZ * M%YY
      CofMt%YY =   M%XX * M%ZZ - M%XZ * M%XZ
      CofMt%YZ = -(M%XX * M%YZ - M%XY * M%XZ)
      CofMt%ZZ =   M%XX * M%YY - M%XY * M%XY
       
      InvertMatS3D = CofMt / DetM
      flops = 35.0
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function InvertMatS3D

#undef __FUNCT__
#define __FUNCT__ "simplexNormal2D"
   !!!
   !!!  
   !!!  simplexNormal2D: Compute the normal to a simplex in 2D
   !!!  
   !!!  (c) 2013 Blaise Bourdin bourdin@lsu.edu
   !!!
   Subroutine simplexNormal2D(Coord,n,ierr)
      Type(Vect2D),Dimension(:),Pointer               :: Coord
      Type(Vect2D),Intent(OUT)                        :: n
      PetscInt      ,Intent(OUT)                      :: ierr

      PetscLogDouble                                  :: flops

      n = (/ Coord(1)%Y-Coord(2)%Y, Coord(2)%X-Coord(1)%X /)
      n = n / norm(n)
      flops = 2
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Subroutine simplexNormal2D

#undef __FUNCT__
#define __FUNCT__ "simplexNormal3D"
   !!!
   !!!  
   !!!  simplexNormal3D: Compute the normal to a simplex in 3D
   !!!  
   !!!  (c) 2013 Blaise Bourdin bourdin@lsu.edu
   !!!
   Subroutine simplexNormal3D(Coord,n,ierr)
      Type(Vect3D),Dimension(:),Pointer               :: Coord
      Type(Vect3D),Intent(OUT)                        :: n
      PetscInt,Intent(OUT)                            :: ierr

      n = (Coord(2)-Coord(1)) .crossP. (Coord(1)-Coord(3))
      n = n / norm(n)
      ierr = 0
   End Subroutine simplexNormal3D

   Function InvertTens4OS2D(T)
      Type(Tens4OS2D),Intent(IN)                  :: T
      Type(Tens4OS2D)                             :: InvertTens4OS2D
      
      PetscInt                                    :: ierr
      PetscReal,Dimension(3,3)                    :: TmpArray
      PetscInt,Dimension(3)                       :: ipiv
      PetscReal,Dimension(3)                      :: work

      !!! We convert T in a matrix using Mandel notations,invert the matrix then write back in a tensor
      TmpArray = T
      Call DGETRF(3,3,TmpArray,3,ipiv,ierr)
      Call DGETRI(3,TmpArray,3,ipiv,work,3,ierr)
      InvertTens4OS2D = TmpArray
   End Function InvertTens4OS2D

   Function InvertTens4OS3D(T)
      Type(Tens4OS3D),Intent(IN)                  :: T
      Type(Tens4OS3D)                             :: InvertTens4OS3D
      
      PetscInt                                    :: ierr
      PetscReal,Dimension(6,6)                    :: TmpArray
      PetscInt,Dimension(6)                       :: ipiv
      PetscReal,Dimension(6)                      :: work

      !!! We convert T in a matrix using Mandel notations,invert the matrix then write back in a tensor
      TmpArray = T
      Call DGETRF(6,6,TmpArray,6,ipiv,ierr)
      Call DGETRI(6,TmpArray,6,ipiv,work,6,ierr)
      InvertTens4OS3D = TmpArray
   End Function InvertTens4OS3D
   
   Subroutine Tens4OS2D2Array4(A,T)
      PetscReal, Dimension(2,2,2,2),Intent(OUT)   :: A
      Type(Tens4OS2D),Intent(IN)                  :: T

      A(1,1,1,1) = T%XXXX
      A(1,1,1,2) = T%XXXY
      A(1,1,2,1) = T%XXXY
      A(1,1,2,2) = T%XXYY
      A(1,2,1,1) = T%XXXY
      A(1,2,1,2) = T%XYXY
      A(1,2,2,1) = T%XYXY
      A(1,2,2,2) = T%YYXY
      A(2,1,1,1) = T%XXXY
      A(2,1,1,2) = T%XYXY
      A(2,1,2,1) = T%XYXY
      A(2,1,2,2) = T%YYXY
      A(2,2,1,1) = T%XXYY
      A(2,2,1,2) = T%YYXY
      A(2,2,2,1) = T%YYXY
      A(2,2,2,2) = T%YYYY
   End Subroutine Tens4OS2D2Array4

   Subroutine Tens4OS3D2Array4(A,T)
      PetscReal, Dimension(3,3,3,3),Intent(OUT)   :: A
      Type(Tens4OS3D),Intent(IN)                  :: T

      A(1,1,1,1) = T%XXXX 
      A(1,1,2,2) = T%XXYY; A(2,2,1,1) = T%XXYY
      A(1,1,3,3) = T%XXZZ; A(3,3,1,1) = T%XXZZ
      A(1,1,2,3) = T%XXYZ; A(1,1,3,2) = T%XXYZ; A(2,3,1,1) = T%XXYZ; A(3,2,1,1) = T%XXYZ;
      A(1,1,1,3) = T%XXXZ; A(1,1,3,1) = T%XXXZ; A(1,3,1,1) = T%XXXZ; A(3,1,1,1) = T%XXXZ;
      A(1,1,1,2) = T%XXXY; A(1,1,2,1) = T%XXXY; A(1,2,1,1) = T%XXXY; A(2,1,1,1) = T%XXXY;

      A(2,2,2,2) = T%YYYY 
      A(2,2,3,3) = T%YYZZ; A(3,3,2,2) = T%YYZZ
      A(2,2,2,3) = T%YYYZ; A(2,2,3,2) = T%YYYZ; A(2,3,2,2) = T%YYYZ; A(3,2,2,2) = T%YYYZ;
      A(2,2,1,3) = T%YYXZ; A(2,2,3,1) = T%YYXZ; A(1,3,2,2) = T%YYXZ; A(3,1,2,2) = T%YYXZ;
      A(2,2,1,2) = T%YYXY; A(2,2,2,1) = T%YYXY; A(1,2,2,2) = T%YYXY; A(2,1,2,2) = T%YYXY;

      A(3,3,3,3) = T%ZZZZ 
      A(3,3,2,3) = T%ZZYZ; A(3,3,3,2) = T%ZZYZ; A(2,3,3,3) = T%ZZYZ; A(3,2,3,3) = T%ZZYZ;
      A(3,3,1,3) = T%ZZXZ; A(3,3,3,1) = T%ZZXZ; A(1,3,3,3) = T%ZZXZ; A(3,1,3,3) = T%ZZXZ;
      A(3,3,1,2) = T%ZZXY; A(3,3,2,1) = T%ZZXY; A(1,2,3,3) = T%ZZXY; A(2,1,3,3) = T%ZZXY;

      A(2,3,2,3) = T%YZYZ; A(2,3,3,2) = T%YZYZ; A(3,2,2,3) = T%YZYZ; A(3,2,3,2) = T%YZYZ;
      A(2,3,1,3) = T%YZXZ; A(2,3,3,1) = T%YZXZ; A(3,2,1,3) = T%YZXZ; A(3,2,3,1) = T%YZXZ;
      A(1,3,2,3) = T%YZXZ; A(1,3,3,2) = T%YZXZ; A(3,1,2,3) = T%YZXZ; A(3,1,3,2) = T%YZXZ;
      A(2,3,1,2) = T%YZXY; A(2,3,2,1) = T%YZXY; A(3,2,1,2) = T%YZXY; A(3,2,2,1) = T%YZXY;
      A(1,2,2,3) = T%YZXY; A(1,2,3,2) = T%YZXY; A(2,3,1,2) = T%YZXY; A(2,3,2,1) = T%YZXY;

      A(1,2,1,2) = T%XZXZ; A(1,2,2,1) = T%XZXZ; A(2,1,1,2) = T%XZXZ; A(2,1,1,2) = T%XZXZ;
      A(1,3,1,2) = T%XZXY; A(1,3,2,1) = T%XZXY; A(3,1,1,2) = T%XZXY; A(3,1,2,1) = T%XZXY;
      A(1,2,1,3) = T%XZXY; A(1,2,3,1) = T%XZXY; A(2,1,1,3) = T%XZXY; A(2,1,3,1) = T%XZXY;
      A(1,2,1,2) = T%XYXY; A(1,2,2,1) = T%XYXY; A(2,1,1,2) = T%XYXY; A(2,1,2,1) = T%XYXY;
   End Subroutine Tens4OS3D2Array4


   Subroutine Array42Tens4OS2D(T,A)
      Type(Tens4OS2D),Intent(OUT)                 :: T
      PetscReal, Dimension(2,2,2,2),Intent(IN)    :: A

      T%XXXX = A(1,1,1,1) 
      T%XXXY = A(1,1,1,2) 
      T%XXYY = A(1,1,2,2) 
      T%XYXY = A(1,2,1,2) 
      T%YYXY = A(2,2,1,2) 
      T%YYYY = A(2,2,2,2) 
   End Subroutine Array42Tens4OS2D

   Subroutine Array42Tens4OS3D(T,A)
      Type(Tens4OS3D),Intent(OUT)                 :: T
      PetscReal, Dimension(3,3,3,3),Intent(IN)    :: A

      T%XXXX = A(1,1,1,1); T%XXYY = A(1,1,2,2); T%XXZZ = A(1,1,3,3); T%XXYZ = A(1,1,2,3); T%XXXZ = A(1,1,1,3); T%XXXY = A(1,1,1,2)
      T%YYYY = A(2,2,2,2); T%YYZZ = A(2,2,3,3); T%YYYZ = A(2,2,2,3); T%YYXZ = A(2,2,1,3); T%YYXY = A(2,2,1,2)
      T%ZZZZ = A(3,3,3,3); T%ZZYZ = A(3,3,2,3); T%ZZXZ = A(3,3,1,3); T%ZZXY = A(3,3,1,2)
      T%YZYZ = A(2,3,2,3); T%YZXZ = A(2,3,1,3); T%YZXY = A(2,3,1,2)
      T%XZXZ = A(1,3,1,3); T%XZXY = A(1,3,1,2)
      T%XYXY = A(1,2,1,2)
   End Subroutine Array42Tens4OS3D

   Function Tens4OS2DTransform(T,M)
      !!! Apply the transformation given by the matrix R to a 4th order tensor
      !!! i.e. C_{ijkl} = R_{ip}.R_{jq}.R_{kr}.R{ls} A_{pqrs}
      Type(Tens4OS2D),Intent(IN)                  :: T
      Type(Mat2D),Intent(IN)                      :: M
      
      Type(Tens4OS2D)                             :: Tens4OS2DTransform
   
      PetscReal,Dimension(2,2,2,2)                :: TT,C
      PetscReal,Dimension(2,2)                    :: MM
      Integer                                     :: i,j,k,l
      Integer                                     :: p,q,r,s
      PetscLogDouble                              :: flops
      PetscErrorCode                              :: ierr
      
      TT = T
      MM = M
      C = 0.0_Kr
      Do i = 1,2
         Do j = 1,2
            Do k = 1,2
               Do l = 1,2
                  Do p = 1,2
                     Do q = 1,2
                        Do r = 1,2
                           Do s = 1,2
                              C(i,j,k,l) = C(i,j,k,l) + MM(i,p) * MM(j,q) * MM(k,r) * MM(l,s) * TT(p,q,r,s)
                           End Do
                        End Do
                     End Do
                  End Do
               End Do
            End Do
         End Do
      End Do
      
      Tens4OS2DTransform = C
      flops = 1280 ! 2**8 * 5
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)

   End Function Tens4OS2DTransform
   
   Function Tens4OS3DTransform(T,M)
      !!! Apply the transformation given by the matrix R to a 4th order tensor
      !!! i.e. C_{ijkl} = R_{ip}.R_{jq}.R_{kr}.R{ls} A_{pqrs}
      Type(Tens4OS3D),Intent(IN)                  :: T
      Type(Mat3D),Intent(IN)                      :: M
      
      Type(Tens4OS3D)                             :: Tens4OS3DTransform
   
      PetscReal,Dimension(3,3,3,3)                :: TT,C
      PetscReal,Dimension(3,3)                    :: MM
      Integer                                     :: i,j,k,l
      Integer                                     :: p,q,r,s
      PetscLogDouble                              :: flops
      PetscErrorCode                              :: ierr

      TT = T
      MM = M
      C = 0.0_Kr
      Do i = 1,3
         Do j = 1,3
            Do k = 1,3
               Do l = 1,3
                  Do p = 1,3
                     Do q = 1,3
                        Do r = 1,3
                           Do s = 1,3
                              C(i,j,k,l) = C(i,j,k,l) + MM(i,p) * MM(j,q) * MM(k,r) * MM(l,s) * TT(p,q,r,s)
                           End Do
                        End Do
                     End Do
                  End Do
               End Do
            End Do
         End Do
      End Do
       
      Tens4OS3DTransform = C
      flops = 32805 ! 3**8 * 5
      Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
   End Function Tens4OS3DTransform

   Function Tens4OS2DSquareRoot(T)
      Type(Tens4OS2D),Intent(IN)                  :: T
      Type(Tens4OS2D)                             :: Tens4OS2DSquareRoot

      Integer,Parameter                           :: n = 3
      Integer                                     :: i,j
      PetscReal,Dimension(n,n)                    :: A,Pt      
      PetscReal,Dimension(n)                      :: lmbda
      PetscReal                                   :: d
      PetscInt                                    :: lwork = 2*n**2+6*n+1
      PetscReal,Dimension(2*n**2+6*n+1)           :: work
      PetscInt                                    :: liwork = 5*n+3
      PetscInt,Dimension(5*n+3)                   :: iwork
      PetscInt                                    :: info

      A = T
      Call DSYEVD('V','L',n,A,n,lmbda,work,lwork,iwork,liwork,info)
      Pt = transpose(A)
      Do i = 1, n
         d = sqrt(lmbda(i))
         If (d < 0.0_Kr) Then
            Write(*,*) 'ERROR in Tens4OSDSquareRoot, negative eigenvalue ', lmbda(i)
         Else
            Do j = 1, n
               Pt(i,j) = A(j,i) * d
            End Do
         End If
      End Do
      Tens4OS2DSquareRoot = matmul(A,Pt)
   End Function Tens4OS2DSquareRoot

   Function Tens4OS3DSquareRoot(T)
      Type(Tens4OS3D),Intent(IN)                  :: T
      Type(Tens4OS3D)                             :: Tens4OS3DSquareRoot

      Integer,Parameter                           :: n = 6
      Integer                                     :: i,j
      PetscReal,Dimension(n,n)                    :: A,Pt      
      PetscReal,Dimension(n)                      :: lmbda
      PetscReal                                   :: d
      PetscInt                                    :: lwork = 2*n**2+6*n+1
      PetscReal,Dimension(2*n**2+6*n+1)           :: work
      PetscInt                                    :: liwork = 5*n+3
      PetscInt,Dimension(5*n+3)                   :: iwork
      PetscInt                                    :: info

      A = T
      Call DSYEVD('V','L',n,A,n,lmbda,work,lwork,iwork,liwork,info)
      Pt = transpose(A)
      Do i = 1, n
         d = sqrt(lmbda(i))
         If (d < 0.0_Kr) Then
            Write(*,*) 'ERROR in Tens4OSDSquareRoot, negative eigenvalue ', lmbda(i)
         Else
            Do j = 1, n
               Pt(i,j) = A(j,i) * d
            End Do
         End If
      End Do
      Tens4OS3DSquareRoot = matmul(A,Pt)
   End Function Tens4OS3DSquareRoot
   
   Subroutine MatS2DSpectralDecomposition(M,ppleValues,ppleDirections)
      Type(MatS2D),Intent(IN)                     :: M
      PetscReal,Dimension(2),Intent(OUT)          :: ppleValues
      Type(MatS2D),Dimension(2),Intent(OUT)       :: ppleDirections
      
      Integer,Parameter                           :: n = 2
      PetscReal,Dimension(n,n)                    :: A,Pt
      Integer                                     :: i
      PetscReal                                   :: d
      PetscInt                                    :: lwork = 2*n**2+6*n+1
      PetscReal,Dimension(2*n**2+6*n+1)           :: work
      PetscInt                                    :: liwork = 5*n+3
      PetscInt,Dimension(5*n+3)                   :: iwork
      PetscInt                                    :: info

      A = M
      Call DSYEVD('V','L',n,A,n,ppleValues,work,lwork,iwork,liwork,info)
      Do i = 1, n
         ppleDirections(i)%XX = A(1,i)**2
         ppleDirections(i)%YY = A(2,i)**2
         ppleDirections(i)%XY = A(1,i) * A(2,i)
      End Do
   End Subroutine MatS2DSpectralDecomposition

   Subroutine MatS3DSpectralDecomposition(M,ppleValues,ppleDirections)
      Type(MatS3D),Intent(IN)                     :: M
      PetscReal,Dimension(3),Intent(OUT)          :: ppleValues
      Type(MatS3D),Dimension(3),Intent(OUT)       :: ppleDirections
      
      Integer,Parameter                           :: n = 3
      PetscReal,Dimension(n,n)                    :: A,Pt      
      Integer                                     :: i
      PetscReal                                   :: d
      PetscInt                                    :: lwork = 2*n**2+6*n+1
      PetscReal,Dimension(2*n**2+6*n+1)           :: work
      PetscInt                                    :: liwork = 5*n+3
      PetscInt,Dimension(5*n+3)                   :: iwork
      PetscInt                                    :: info

      A = M
      Call DSYEVD('V','L',n,A,n,ppleValues,work,lwork,iwork,liwork,info)
      Do i = 1, n
         ppleDirections(i)%XX = A(1,i)**2
         ppleDirections(i)%YY = A(2,i)**2
         ppleDirections(i)%ZZ = A(3,i)**2
         ppleDirections(i)%YZ = A(2,i) * A(3,i)
         ppleDirections(i)%XZ = A(1,i) * A(3,i)
         ppleDirections(i)%XY = A(1,i) * A(2,i)
      End Do
   End Subroutine MatS3DSpectralDecomposition



   Subroutine MatS3DEigenVectorValues(M,MatProj,MatDiag)
      Type(MatS3D),Intent(IN)                     :: M
      PetscReal,Dimension(3)                      :: ppleValues
      Type(Mat3D),Intent(OUT)                     :: MatProj
      Type(MatS3D),Intent(OUT)                    :: MatDiag
      
      Integer,Parameter                           :: n = 3
      PetscReal,Dimension(n,n)                    :: A,Pt      
      Integer                                     :: i
      PetscReal                                   :: d
      PetscInt                                    :: lwork = 2*n**2+6*n+1
      PetscReal,Dimension(2*n**2+6*n+1)           :: work
      PetscInt                                    :: liwork = 5*n+3
      PetscInt,Dimension(5*n+3)                   :: iwork
      PetscInt                                    :: info

      A = M
      Call DSYEVD('V','L',n,A,n,ppleValues,work,lwork,iwork,liwork,info)

      MatDiag=0.0_Kr
      MatDiag%XX=ppleValues(1)
      MatDiag%YY=ppleValues(2)
      MatDiag%ZZ=ppleValues(3)

      MatProj%XX = A(1,1)
      MatProj%YX = A(2,1)
      MatProj%ZX = A(3,1)
      MatProj%XY = A(1,2)
      MatProj%YY = A(2,2)
      MatProj%ZY = A(3,2)
      MatProj%XZ = A(1,3)
      MatProj%YZ = A(2,3)
      MatProj%ZZ = A(3,3)
   End Subroutine MatS3DEigenVectorValues


   Subroutine MatS2DEigenVectorValues(M,MatProj,MatDiag)
      Type(MatS2D),Intent(IN)                     :: M
      PetscReal,Dimension(2)                      :: ppleValues
      Type(Mat2D),Intent(OUT)                     :: MatProj
      Type(MatS2D),Intent(OUT)                    :: MatDiag
      
      Integer,Parameter                           :: n = 2
      PetscReal,Dimension(n,n)                    :: A,Pt      
      Integer                                     :: i
      PetscReal                                   :: d
      PetscInt                                    :: lwork = 2*n**2+6*n+1
      PetscReal,Dimension(2*n**2+6*n+1)           :: work
      PetscInt                                    :: liwork = 5*n+3
      PetscInt,Dimension(5*n+3)                   :: iwork
      PetscInt                                    :: info
      PetscErrorCode                              :: ierr

      A = M
      Call DSYEVD('V','L',n,A,n,ppleValues,work,lwork,iwork,liwork,info)
      If (info /= 0) Then
         Write(*,*) 'DSYEVD failed with info=',info
         Write(*,*) 'A: ', A
         Write(*,*) 'ppleValues: ',ppleValues
         SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_LIB,"DSYEVD failed: "//__FUNCT__,ierr)
      End If

      MatDiag=0.0_Kr
      MatDiag%XX=ppleValues(1)
      MatDiag%YY=ppleValues(2)

      MatProj%XX = A(1,1)
      MatProj%XY = A(1,2)
      MatProj%YX = A(2,1)
      MatProj%YY = A(2,2)
   End Subroutine MatS2DEigenVectorValues

#undef __FUNCT__
#define __FUNCT__ "Mat2DMoment"
!!!
!!!  
!!!  Mat2DMoment: k-th moment of a 2x2 matrix, i.e tr(A^k)
!!!  
!!!  (c) 2019 Blaise Bourdin bourdin@lsu.edu
!!!

   Function Mat2DMoment(k,A)
      Integer,Intent(IN)                          :: k
      Type(Mat2D),Intent(IN)                      :: A
      PetscReal                                   :: Mat2DMoment

      Integer                                     :: i
      Type(Mat2D)                                 :: Ak
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      Select Case(k)
      Case(1)
         Mat2DMoment = trace(A)
         flops = 0.
      Case(2)
         Mat2DMoment = (A%XX**2 + 2*A%XY*A%YX + A%YY**2) / 2.0_Kr
         flops = 7.
         Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
      Case default
         Ak = A
         do i = 1, k-1
            Ak = Ak*A
         end do 
         Mat2DMoment = trace(Ak)/k
         flops = 1.
      End Select
   End Function Mat2DMoment

#undef __FUNCT__
#define __FUNCT__ "MatS2DMoment"
!!!
!!!  
!!!  MatS2DMoment: k-th moment of a 2x2 symmetric matrix, i.e tr(A^k)
!!!  
!!!  (c) 2019 Blaise Bourdin bourdin@lsu.edu
!!!

   Function MatS2DMoment(k,A)
      Integer,Intent(IN)                          :: k
      Type(MatS2D),Intent(IN)                     :: A
      PetscReal                                   :: MatS2DMoment

      Integer                                     :: i
      Type(MatS2D)                                :: Ak
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      Select Case(k)
      Case(1)
         MatS2DMoment = trace(A)
         flops = 0.
      Case(2)
         MatS2DMoment = (A%XX**2 + 2*A%XY**2 + A%YY**2) / 2.0_Kr
         flops = 7.
         Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
      Case default
         Ak = A
         do i = 1, k-1
            Ak = Ak*A
         end do 
         MatS2DMoment = trace(Ak) / k
         flops = 1.
      End Select
   End Function MatS2DMoment

#undef __FUNCT__
#define __FUNCT__ "Mat3DMoment"
!!!
!!!  
!!!  Mat3DMoment: k-th moment of a 3x3 matrix, i.e tr(A^k)
!!!  
!!!  (c) 2019 Blaise Bourdin bourdin@lsu.edu
!!!

   Function Mat3DMoment(k,A)
      Integer,Intent(IN)                          :: k
      Type(Mat3D),Intent(IN)                      :: A
      PetscReal                                   :: Mat3DMoment

      Integer                                     :: i
      Type(Mat3D)                                 :: Ak
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      Select Case(k)
      Case(1)
         Mat3DMoment = trace(A)
         flops = 0.
      Case(2)
         Mat3DMoment = (A%XX**2 + A%YY**2 + A%ZZ**2 + 2*A%XY*A%YX + 2*A%XZ*A%ZX + 2*A%YZ*A%ZY) / 2.0_Kr
         flops = 15.0
         Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
      Case default
         Ak = A
         do i = 1, k-1
            Ak = Ak*A
         end do 
         Mat3DMoment = trace(Ak) / k
         flops = 1.
      End Select
   End Function Mat3DMoment

#undef __FUNCT__
#define __FUNCT__ "MatS3DMoment"
!!!
!!!  
!!!  MatS3DMoment: k-th moment of a 3x3 symmetric matrix, i.e tr(A^k)
!!!  
!!!  (c) 2019 Blaise Bourdin bourdin@lsu.edu
!!!

   Function MatS3DMoment(k,A)
      Integer,Intent(IN)                          :: k
      Type(MatS3D),Intent(IN)                     :: A
      PetscReal                                   :: MatS3DMoment

      Integer                                     :: i
      Type(MatS3D)                                :: Ak
      PetscLogDouble                              :: flops
      PetscInt                                    :: ierr
      Select Case(k)
      Case(1)
         MatS3DMoment = trace(A)
         flops = 0.
      Case(2)
         MatS3DMoment = (A%XX**2 + A%YY**2 + A%ZZ**2 + 2*A%YZ**2 + 2*A%XZ**2 + 2*A%XY**2) / 2.0_Kr
         flops = 15.
         Call PetscLogFlops(flops,ierr);CHKERRQ(ierr)
      Case default
         Ak = A
         do i = 1, k-1
            Ak = Ak*A
         end do 
         MatS3DMoment = trace(Ak) / k
         flops = 1.
      End Select
   End Function MatS3DMoment

End Module m_MEF90_LinAlg
